'  AD-A129  461  THE  NUMERICAL  SYNTHESIS  AND  INVERSION  OF  ACOUSTIC  1/t 

FIELDS  USING  THE  HANKEL . . (U)  MASSACHUSETTS  INST  OF  TECH 
CAMBRIDGE  RESEARCH  LAB  OF  ELECTRON..  D  R  MOOK  JAN  83 
UNCLASSIFIED  TR-497  N00014-81-K-0742  F/G  20/1  NL 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  Of  STANDARDS-1963-A 


LE  COPY 


"The  Numerical  Synthesis  and  Inversion  of 
Acoustic  Fields  Using  the  Hankel  Transform 
Q  with  Application  to  the  Estimation  of  the 

•^C  Plane  Wave  Reflection  Coefficient  of  the 

Ocean  Bottom” 

Douglas  R.  Mepk 

Technical  Report  497  -  January  1983 
M.I.T.,  R.L.E. 


DTIC 


This  work  has  bean  supported  in  part  by  the  Advanced  Research  Projects  Agency 
monitored  by  ONR  under  Contract  N00014-81-K-0742  NR-049-506  and  in  part  by  the 
National  Science  Foundation  under  Grant  ECS80-07102. 


^Public  leiww  and 
distribution  j3  uniimit^  '  Us 


SECURITY  CLASSIFICATION  OF  THIS  A  AGE  (Whit  Dole  Cnttctd) 


|  REPORT  DOCUMENTATION  PAGE 

SEAS  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

r . . 

4.  TITLE  (md  Submit) 

The  Numerical  Synthesis  and  Inversion  of 
Acoustic  Fields  Using  the  Hankel  Transform 
with  Application  to  the  Estimation  of  the 
Plane  Wave  Reflection  Coefficient  of  the  O' 

S.  TYPE  or  REPORT  4  PERIOD  COVERED 

Technical  Report 

S.  PERFORMING  ORO.  REPORT  NUMBER 

:ean  Bottom 

7.  AuTMORf*) 

Douglas  R.  Mook 

I.  CONTRACT  OR  GRANT  HUMBERT*.) 

N00014-81-K-0742 

S.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Research  Laboratory  of  Electronics 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139 

10.  program  element,  project.  TASK 
AREA  •  WORK  UNIT  NUMBERS 

NR- 049-506 

II.  CONTROLLING  OFFICE  NAME  AND  AOONESS 

Advanced  Research  Projects  Agency 

1400  Wilson  Boulevard 

Arlington.  VA  22217  . , 

It.  REPORT  DATE 

January  1983 

IS.  number  of  RAGES 

226 

14.  MONITORING  agency  NAME  4  AOORESS flldllltrmt  from  Controlling  OUlet) 

Office  of  Naval  Research 

Mathematical  and  Information  Sciences 
Division,  800  North  Quincy  Street 

Arlington,  VA  22217 

is.  security  class,  toi  mi*  report) 

Unclassified 

IS*.  DECLASSIFICATION  downgrading 

schedule 

14.  DISTRIBUTION  STATEMENT  (ol  tbit  Report) 

approved  for  public  release;  distribution  unlimited 

17.  DISTRIBUTION  STATEMENT  (of  tbt  tbatrtet  tnttrtd  In  Block  20,  II  d‘  ■  »r#ni  bom  Report) 

IS.  SUPPLEMENTARY  NOTES 

Hankel  transform  plane  wave  reflection  coefficient 
Sommerfeld  integral  depth-dependent  Green's  function 
trapped  waves  acoustic  CW  point  source 

20.  ABSTRACT  fConitou*  an  rtvtrtt  tldt  II  ntctttory  md  Idtntllr  by  block  mrnifcar) 

see  other  side 

oo  ,:2r„  W73 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OP  THIS  PAGE  3»SI  Dote  Entered) 


UNCLASSIFIED _ 

SSCUWITV  CLASSIFICATION  or  THIS  VAgKfWhn  Data  hum) 


ABSTRACT 


The  plane  wave  reflection  coefficient  is  an  important  geometry  independent  means  of  speci¬ 
fying  the  acoustic  response  of  a  horizontally  stratified  ocean  bottom.  It  is  an  integral  step  in  file 
inversion  of  acoustic  field  measurements  to  obtain  parameters  of  the  bottom  and  it  is  used  to 
characterize  an  environment  for  purposes  of  acoustic  imaging.  This  diesis  studies  both  the  gen¬ 
eration  of  synthetic  pressure  fields  through  the  plane  wave  reflection  coefficient  and  the  inversion 
of  measured  pressure  fields  to  estimate  die  plane  wave  reflection  coefficient.  These  are  related 
through  the  Sommexfeld  integral  which  is  in  the  form  of  a  Hankcl  transform.  The  Hankel 
transform  is  extensively  studied  in  this  diesis  and  both  theoretical  properties  and  numerical 
implementations  are  considered.  These  results  have  broad  applications.  When  we  apply  diem  to 
the  generation  of  synthetic  data,  we  obtain  hybrid  numerical-analytical  algorithms  which  pro¬ 
vide  extremely  accurate  synthetic  fields  without  sacrifiring  computational  speed.  These  algorithms 
can  accurately  incorporate  die  effects  of  trapped  modes  guided  by  slow  speed  layers  in  die  bot¬ 
tom.  We  also  apply  these  tools  to  study  the  inversion  of  measured  pressure  field  data  for  die 
plane  wave  reflection  coefficient.  We  address  practical  issues  associated  with  the  inversion  pro¬ 
cedure  including  removal  of  the  source  field,  sampling,  field  measurements  over  a  finite  range, 
and  uncontrolled  variations  in  source-height  A  phase  unwrapping  and  associated  interpolation 
scheme  is  developed  to  handle  improperly  spaced  data. 

A  preliminary  inversion  of  real  pressure  field  data  is  performed.  In  parallel,  an  inversion 
of  a  synthetically  generated  field  for  similar  bottom  parameters  is  also  performed  and  the  results 
of  processing  the  real  and  synthetic  data  are  compared.  The  estimate  for  the  depth  dependent 
Green’s  function  obtained  from  the  real  data  shares  many  features  with  die  depth  dependent 
Green’s,  function  estimated  from  the  synthetic  data,  suggesting  that  the  total  inversion  to  obtain 
the  plane  wave  reflection  coefficient  will  soon  be  possible.  Errors  in  die  present  estimate  of  the 
plane  wave  reflection  coefficient  are  associated  with  uncontrolled  source-height  variations  during 
the  acquisition  of  data. 
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The  plane  wave  reflection  coefficient  is  an  important  geometry  independent  means  of  speci¬ 
fying  the  acoustic  response  of  a  horizontally  stratified  ocean  bottom.  It  is  an  integral  step  in  the 
inversion  of  acoustic  field  measurements  to  obtain  parameters  of  the  bottom  and  it  is  used  to 
characterize  an  environment  for  purposes  of  acoustic  imaging.  This  thesis  studies  both  the  gen¬ 
eration  of  synthetic  pressure  fields  through  the  plane  wave  reflection  coefficient  and  the  inversion 
of  measured  pressure  fields  to  estimate  die  plane  wave  reflection  coefficient. .  These  are  related 
through  the  Sommerfdd  integral  which  is  in  the  form  of  a  Hankel  transform.  The  Hankel 
transform  is  extensively  studied  in  this  thesis  and  both  theoretical  properties  and  numerical 
implementations  are  considered.  These  remits  have  broad  applications.  When  we  apply  them  to 
the  generation  of  synthetic  data,  we  obtain  hybrid  numerical-analytical  algorithms  which  pro¬ 
vide  extremely  accurate  synthetic  fields  without  sacrifinng  computational  speed.  These  algorithms 
can  accurately  incorporate  die  effects  of  trapped  modes  guided  by  slow  speed  layers  in  die  bot¬ 
tom.  We  also  apply  these  tools  to  study  die  inversion  of  measured  presrare  field  data  for  die 
plane  wave  reflection  coefficient.  We  address  practical  isroes  amodated  with  the  inversion  pro¬ 
cedure  including  removal  of  die  source  field,  sampling,  field  measurements  over  a  Gnu*  range, 
and  uncontrolled  variations  in  source-height.  A  phase  unwrapping  and  associated  interpolation 
scheme  is  developed  to  handle  improperly  spaced  data. 

A  preliminary  inverrion  of  real  pressure  field  data  is  performed.  In  parallel,  an  inversion 
of  a  synthetically  generated  field  for  similar  bottom  parameters  is  also  berforfaed  and  the  results 
of  processing  die  real  and  synthetic  data  are  compared.  The  estimate!  fohfhe  depth  dependent 
Green's  function  obtained  from  the  real  data  shares  many  features  with  me  depth  dependent 
Green's  function  estimated  from  the  synthetic  data,  suggesting  that  die  total  inversion  to  obtain 
the  plane  wave  reflection  coefficient  will  soon  be  possible.  Errors  in  the  presentsestimate  of  the 
plane  wave  reflection  coefficient  are  associated  with  uncontrolled  source-height  vanathms  during 
the  acquisition  of  data.  \ 
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CHAPTER  I: 

INTRODUCTION 

1.1)  Overview 

The  plane  wave  reflection  coefficient  is  an  important  geometry  independent  means  of  speci¬ 
fying  die  acoustic  response  of  a  horizontally  stratified  ocean  bottom.  It  is  used  both  to  cnlailat«» 
fields  in  propagations  models  and  as  an  input  to  a  variety  of  inverse  techniques  which  seek  to 
determine  bottom  parameters  [1,2,3  ].  In  this  diesis  we  will  study  both  die  generation  of  syn¬ 
thetic  pressure  fields  through  the  plane  wave  reflection  coefficient  and  the  inversion  of  measured 
pressure  fields  to  estimate  the  plane  wave  reflection  coefficient.  We  will  consider  only  die  fields 
associated  with  a  CW  point  source  in  the  deep  ocean  over  a  horizontally  stratified  bottom  and 
will  not  allow  the  bottom  to  support  shear  waves.  The  results  of  this  thesis,  however,  are  appli¬ 
cable  to  a  wide  dass  of  wave  problems  and  can  be  generalized  to  permit  the  source  to  be  within 
any  isovelocity  layer  and  the  introduction  of  shear.  Further,  it  is  our  hope  that  the  tools 
developed  in  the  course  of  this  work  will  find  applications  in  many  areas. 

The  foundation  for  our  studies  of  the  forward  and  inverse  problems  is  die  Hankel 
transform,  [4, 5, 6  ]  which  arises  in  these  contexts  because  the  Sommerfeld  integral  which  relates 
the  plane  wave  reflection  coefficient  to  the  reflected  field  is  in  that  form.  We  wifi  derive  general 
properties  of  the  Hankel  transform  to  guide  the  work  in  these  areas.  We  will  also  study  and 
develop  numerical  implementations  to  permit  computer  processing.  A  fast,  accurate  numerical 
Hankel  transform  algorithm  is  developed  and  illustrated. 

The  Hankel  transform  results  allow  us  to  isolate  significant  sources  of  degradation  in 
numerically  generated  synthetic  fields.  To  remove  these  sources  we  develop  a  hybrid  analytical- 
numerical  procedure  which  is  significantly  more  accurate  without  sacrificing  computational 
speed.  This  hybrid  algorithm  performs  some  of  the  calculations  analytically  while  keeping  the 
numerical  computation  in  the  form  of  a  Hankel  transform. 

Through  another  hybrid  technique  we  incorporate  the  effects  of  trapped  modes  diet  may  be 
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guided  by  low  speed  layers  in  the  bottom.  The  reaults  are  accurate  both  in  the  near  and  far 
fields,  in  contrast  with  modal  methods.  The  method  developed  for  handling  these  modes  can 
also  be  used  to  control  other  complications  associated  with  poles  in  the  plane  wave  reflection 
coefficient  such  as  those  that  would  be  introduced  by  allowing  for  shear  wave  propagation. 

In  Chapter  V  we  begin  to  study  the  inversion  of  pressure  field  data  to  obtain  an 
for  the  plane  wave  reflection  coefficient  We  draw  upon  our  previous  results  to  consider  a 
recently  proposed  method  for  performing  this  inversion  by  Frisk,  Oppenheim  and  Martinez  [7  ]. 
Frisk,  Oppenheim  and  Martinez  have  proposed  that  the  Sommerfeld  integral  be  inverted 
directly,  without  recourse  to  the  specular  approximation  used  by  previous  methods.  [7  ]  With 
such  a  direct  inversion,  estimates  would  no  longer  be  confined  to  real  angles  of  incidence  and 
regions  where  the  specular  approximation  is  valid.  [8, 9  ]  Such  a  direct  inversion  has  been  made 
possible  recently  by  coherent  measurements  of  the  reflected  presmre  field  resulting  from  a  point 
source  over  a  horizontally  stratified  bottom.  [10  ]  In  this  chapter  we  study  several  practical  issues 
associated  with  this  proposed  direct  inversion.  We  consider  directly  the  issues  of  source  field  sub¬ 
traction,  sampling,  windowing  of  the  pressure  field,  and  uncontrolled  variations  in  source  height 
The  issue  of  source-field  subtraction  arises  because  the  plane  wave  reflection  coefficient  is 
directly  related  to  the  reflected  field  and  not  the  total  field,  which  is  measured.  Under  the  issue 
of  sampling  we  study  both  the  sampling  rate  required  to  obtain  a  valid  inversion  and  the  effects 
of  improperly  spaced  data.  To  handle  improperly  spaced  data,  an  interpolation  procedure  is 
developed  which  is  based  on  a  new  phase  unwrapping  procedure.  In  the  discussion  of  window¬ 
ing  we  determine  the  range  to  which  field  measurements  must  be  taken  in  order  to  obtain  a  valid 
inversion,  hi  the  section  on  source-height  variation  we  characterize  the  degradation  that  results 
from  variation  in  the  source-height  and  study  canonical  variations. 

Having  considered  many  of  the  isnes  affecting  the  direct  inversion  of  pressure  field  data  to 
obtain  the  plane  wave  reflection  coefficient  we  now  actually  perform  a  preliminary  inverson  of 
real  data.  In  parallel  we  invert  synthetic  data  that  we  have  generated  unng  bottom  parameters 


comparable  to  thorn  we  believe  associated  with  the  real  data.  We  draw  upon  the  previous 
developments  to  interpret  the  results. 

In  the  remaining  portion  of  this  diapter  we  briefly  develop  the  acoustic  framework  upon 
which  this  thesis  rests.  We  also  describe  the  experimental  paradigm  by  which  the  real  data  was 
gathered. 

L2)  Plant  Waves  and  a  Horizontally  Stratified  Environment 

A  horizontally  stratified  environment  is  one  for  which  file  material  parameters  vary  only 
vertically.  Such  a  simple  environment  makes  posable  an  in  depth  study  without  the  complications 
that  a  more  varied  environment  would  introduce.  The  insights  gained  from  studying  this  simple 
environment  can  provide  an  understanding  of  more  complicated  envonments.  Also,  for  many 
conditions,  such  as  are  found  in  the  region  of  an  abyssal  plane  in  the  deep  ocean,  the  model  is 
itself  sufficient. 

Figure  1.2.1  shows  an  isovelodty  water  layer  over  a  horizontally  stratified  bottom.  Within 
the  water  a  single  plane  wave  has  the  form  :* 

+4,*)^  —  lut 

For  this  wave  to  be  a  solution  to  the  wave  equation  within  the  water  (which  has  sound  speed  c  ): 

v*_7'?r  w 

the  wave  numbers  (kx ,  kf ,  and  kt)  must  satisfy: 

-  4  -  0  (3) 

c 

We  define  the  vector  k  m  kxtx  +kyty  +kxix  and  the  scalar,  k  ■  k  will  point  in  die  direc¬ 
tion  that  the  plane  wave  propagates.  In  terms  of  this  vector,  the  requirement  (3)  can  be  written 

11*11  -  Vk*+k*+k*  -  *  (4) 

1)  Throughout  this  that*  w»  win  ai pprm  the  *  ~l*  Ha*  daptadana  bacaua  It  uparstw  out  froa  an 
tta  fluid  axpraaUoas. 
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Figure  1.2.1  Incident  plane  wave  geometry 


-  14  - 


If  (4)  is  satisfied  then  (1)  is  formally  a  solution  to  the  wave  equation  even  if  £  is  complex.  When 
one  or  more  of  the  components  of  k  are  imaginary  the  plane  wave  is  called  evanescent  and  will 
vary  exponentially  in  the  direction  of  the  imaginary  component^). 

Considerations  of  symmetry  guarantee  that  when  a  plane  wave  strikes  a  horizontally  strati' 
fied  bottom  the  resulting  wave  will  also  be  planar.  Fc?  the  purposes  of  field  calculations,  plane 
waves  are  eigenfunctions  of  horizontally  stratified  systems.  An  incident  plane  wave  given  by  : 

(5) 

will  generate  the  reflected  wave 

p*ei(k,x+kfJ-k,z)  (g) 

The  change  in  sign  indicates  that  the  reflected  wave  is  returning  in  the  z  direction.  The  ampli¬ 
tude  change  defines  the  plane  wave  reflection  coefficient,  which  may  be  a  function  of  k .  Since 
for  fixed  u>  only  two  of  the  three  components  of  k  can  be  specified  independently  we  write  the 
plane  wave  reflection  coefficient  explicitly  as  a  function  of  only  two: 

P> 

Our  plane  wave  reflection  coefficient  is  defined  for  a  single  frequency  only.  This  implies 
that  the  incident  signal  has  been  present  for  all  time.  The  returning  wave,  PRe^k*x  +k,}  k‘*\  is 
influenced  by  the  bottom  at  all  depths.  This  is  in  contrast  to  the  occasional  usage  of  the  term  for 
which  only  the  surface  contribution  to  a  pulsed  input  is  considered. 

1.3)  The  Weyl  Integral 

Because  die  wave  propagation  we  consider  is  linear  we  can  determine  the  response  to  a 
more  complicated  incident  field  by  considering  that  field  as  a  superpoation  of  plane  waves.  [9  ] 

To  calculate  the  reflected  responses  to  a  point  source  shown  in  Figure  1.3.1  we  first  express 

the  fold  of  the  point  source  at  x  •  zq  as  a  superposition  of  plane  waves.  We  write:1 
1)  Wewifl  use  ftp  to  denote  fee  wave  number  in  the  wear. 
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Figure  1.3.1  Point  source  geometry 


The  field  associated  with  a  point  source  over  a  horizontally  stratified  bottom  is  symmetric 
about  the  z  axis  of  Figure  1.3.1.  We  can  exploit  this  symmetry  to  reduce  the  two  dimensional 
Weyi  integral  (1.3.5)  to  a  one  dimensional  integral. 

The  symmetry  of  the  problem  guarantees  that  all  die  field  variables  in  space  show  a  cylindr¬ 
ical  symmetry.  Because  the  two  Himmamul  Fourier  transform  of  a  circularly  symmetric  function 
win  also  be  drculady  symmetric,  the  Fourier  domain  wiQ  also  display  a  cylindrical  symmetry. 
We  define 


r 2  m  jr2+y2 


k*  - 

With  these  definitions  we  write  (1.3.5)  in  cylindrical  coordinates  as: 

=  T7//77-? - 

**0  0  vkQ—kf 

With Tp/kg+ky)  m  rc (kx,kj).  Performing  the  8  integration  and  using  [14  ] 


lAx)  - 

zir  0 


Equation  (2)  becomes 


(1) 


(2) 


(3) 


F*(r,z)  »  r(k,)e,V**^H't+')J0(k,r)krdkr 

This  is  the  Sommerfdd  integral.  [13  ] 

The  Sommerfeld  integral  has  the  form 


(4) 


P*(r,t)  «  jG{krtz^0y^rkr)krdkr  (5) 

0 

where 

We  wfil  refer  to  G(krrx  jq)  as  die  depth  dependent  Green’s  function  or  the  Green’s  function. 
The  integral  transform  (5)  is  die  Hankei  transform.  [4, 5, 6  ] 


} 

\ 
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Equation  (4)  represents  the  reflected  field  as  a  superposition  of  cylindrical  waves  of  the 

form  J  oCkf  r  )e  \  eadi  of  which  can  be  considered  to  be  a  superposition  of  plane 

waves  all  striirin|  the  surface  with  the  same  horizontal  wave  number,  kr.  Because  T(hr)  is 
r dated  so  directly  to  die  plane  wave  reflection  coefficient,  T(\/k2+k2)  -  rc (kxtkj),  we  wil 
refer  to  it  as  the  plane  wave  reflection  coefficient 


I.S)  Obtaining  the  Reflection  Coefficient 

Presently  most  techniques  for  determining  die  plane  wave  reflection  coefficient  as  a  func¬ 
tion  of  horizontal  wave  number,  r(f,),  from  die  reflected  pressure  fidd,  Pg  (r),  do  not  concen¬ 
trate  on  inverting  Equation  (1.4.4).  Instead  they  consider  the  stationary  phase  approximation  to 
that  Equation  given  by:1  [9  ] 


r(*o)etto* 

Pjt(')  * - ^ - 

where  k0  is  defined  to  be  the  water  wave  number,  — ,  and  R2  m  r2+(z0+z )2 

c 

Equation  (L4.4)  is  inverted  to  provide 


(1) 


r<X>= 


“o* 


PK(r) 


(2) 


S  Equation  (2)  is  used  to  estimate  T(-— -)  then  |r(—~ )  |  can  be  from  the  mag- 

K  K 


nitude  of  Pr(r)  alone  through: 


u«i  moi  ro 

This  1m ct  together  with  the  simplicity  of  Equation  (2)  account  for  its  widespread  use. 


I)  Mew  tin— iyAsplsaswavs  raflsstlns  ssdBsAat  m  s  tenet Isa  at  wfe  Is  ralswd  Aenagh  As  tgem- 
he  sppraAnsdsa  is  hAnM  pnmne  Ml  Ws  eeet  >  ss  s  fsanrtna  of  hortasaal  wm  a—bsrto  be 
tmetemm  et*  As  ns  at  As  ubl  If  ws  tfssow  r#(9)  «  As  plsss  «m  reBseUea  codMsat  ss  s  (u» 

Use  sf  sagh.  As  tee  AntsMses  an  rstissd  Aroagfc  T«(t)  ■  r(ftgsin(9)).  i(j  is  As  horisoatal 


A  As 


Unfortunately  the  stationary  phase  approximation  upon  which  (4)  is  based  is  appropriate  only 
for  distances,  R ,  large  compared  to  a  wavelength  and  only  for  specular  angles  less  than  criti¬ 
cal.^, 9]  It  completely  ignores  near  field  effects  associated  with  T(kr)  for  kr  >ko.  For  applica¬ 
tions  that  consider  near  field  effects  or  the  character  of  the  pressure  field  close  to  or  greater  than 
the  critical  angle,  a  more  exact  invernon  of  Equation  (1.4.4)  is  required. 

For  such  applications  Frisk,  Oppenheim,  and  Martinez  [7]  have  proposed  that  both  the 
magnitude  and  phase  of  the  reflected  pressure  field  be  measured  and  that  the  Sommerfeld 
integral  of  Equation  (1.4.4)  be  inverted  directly,  without  recourse  to  the  stationary  phase 
approximation.  The  Sommerfeld  integral  is  in  the  form  of  a  Hank  el  transform.  Sm«  the 
Hank  el  transform  is  its  own  inverse  [5],  Equation  (1.4.4)  can  be  inverted  and  solved  for  the 
plane  wave  reflection  coefficient  in  terms  of  Pg(r),  the  reflected  response  to  a  paint  source. 
This  gives: 

F(*,)  «  -iVki-kr2e~iV^  ‘  +Z°  fPjt(r )J0(krr )rdr  (5) 

0 


1.6)  Experimental  Model 

In  this  thesis  we  wifi  perform  a  preliminary  inversion  of  measured  pressure  field  data 
according  to  Equation  (1.5.5),  as  suggested  by  Frisk,  Oppenheim,  and  Martinez.  [7]  The  data 
we  analyze  was  taken  by  O.  Frisk,  J.  Doutt,  and  E.  Hays  in  1981.  In  this  section  we  present  the 
details  of  the  experimental  configuration  they  used.  A  experimental  configuration  has 

been  described  in  the  literature.  [10] 

Figure  L6.1  tiiows  the  experimental  configuration  used  by  Frisk,  Doutt,  and  Hays.  As 
feown,  two  receivers  were  moored  1.17  and  54.55  meters  from  the  bottom  of  the  ocean  on  an 
abyssal  plain  under  1900  meters  of  ocean.  The  source  was  attached  by  cable  to  a  drip  on  fee  sur¬ 
face,  and  drifted  slowly  away  at  a  height  off  the  bottom  of  approximately  135  meters.  Every  12 
seconds  the  source  emitted  a  4  second  220  Hertz  tone  which  fee  receivers  recorded  after  quadra- 
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tuze  demodulation  and  low  pan  filtering,  fit  this  way  the  complex  amplitude  of  the  field  as  a 
function  of  range  was  recorded.  The  source  strength  was  177  dB  re  1  fiPa  at  1  meter. 

Recording  by  the  receivers  was  initiated  every  12  seconds  upon  recognition  of  an  11  kHz 
trigger  pulse  sent  from  a  pinger  mounted  on  die  source,  and  continued  for  6  seconds.  During  this 
time  the  output  of  the  receivers  was  quadrature  demodulated,  low  pan  filtered  to  2  Hz,  digitized 
by  a  12  bit  A/D  converter  at  a  5  Hz  rate  and  recorded  on  cassette  tape.  A  schematic  of  the  prel¬ 
iminary  data  processing  in  die  receiver  system  is  shown  in  Figure  1.6.2. 

The  ship  drifted  at  a  speed  of  about  1/2  kn  allowing  one  sample  of  the  field  every  half 
wavelength.  The  clocks  in  the  source  and  receiver  were  synchronized  and  had  a  stability  of  about 
one  part  in  109  per  day.  The  11  kHz  emission  times  at  the  source  and  the  arrival  time  at  the 
receivers  were  used  to  determine  the  slant  range  between  die  source  and  the  receivers.  As  part  of 
die  processing  it  was  necessary  to  estimate  the  source  height  and  convert  from  riant  ranged  to 
horizontal  range. 

Frisk,  Doutt  and  Hays  determined  that  the  signal  was  in  a  steady  state  condition  by  the  4th 
data  sample. 

1.7)  Summary 

In  this  thesis  we  consider  the  generation  of  synthetic  pressure  fields  through  the  evaluation 
of  the  Sommerfeld  integral.  This  integral  is  in  the  form  of  the  Hankel  transform  of  the  depth- 
dependent  Green's  function.  We  also  consider  the  inversion  of  a  measured  pressure  field  to  esti¬ 
mate  the  depth-dependent  Green's  function  and  from  that  the  plane  wave  reflection  coefficient 
The  foundation  of  both  these  procedures  is  die  Hankel  transform.  In  die  next  chapter  we  will 
both  catalogue  and  develop  the  properties  of  the  Hankel  transform  that  will  provide  the  founda¬ 
tion  for  the  work  of  this  thesis. 
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Fignre  1.6.2  Schematic  of  preliminary  processing  by  data  acquisition  system 
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CHAPTER  n: 

THE  HANKEL  TRANSFORM 

n.l)  Overview 

The  relation  of  the  Hankel  transform  to  the  two  dimensional  Fourier  transform  of  a  circu¬ 
larly  symmetric  function  makes  it  as  important  a  tool  for  problems  cast  in  cylindrical  coordinate 
systems  as  die  Fourier  transform  for  proplems  in  cartesian  systems.  Applications  can  be  found  in 
such  diverse  fields  as  astronomy,  electrodynamics,  electrostatics,  oceanography,  physics,  and 
seismology.  Because  it  relates  the  pressure  field  associated  with  a  point  source  in  a  horizontally 
stratified  medium  to  the  plane  wave  reflection  coefficient,  it  forms  the  foundation  of  this  thesis. 
In  this  chapter  we  explore  the  properties  of  the  Hankel  transform. 

We  begin  by  presenting  the  most  common  definitions  of  the  Hankel  transform  in  Section 
II.2.  We  show  how  the  Hankel  transform  arises  from  the  two  dimensional  Fourier  transform  of  a 
circularly  symmetric  function  in  Section  II. 3.  To  relate  the  Hankel  transform  to  the  more  fami¬ 
liar  one  dimensional  Fourier  transform,  in  Section  Ct.4  we  present  its  asymptotic  form.  In  Section 
n.5  we  complete  our  presentation  of  available  properties  with  a  summary  of  important  results 
available  in  the  literature. 

The  remainder  of  this  chapter  is  devoted  to  results  previously  unavailable.  We  derive  these 
results  to  provide  the  foundation  for  our  work  later  in  this  thesis.  Section  II. 6  examines  window¬ 
ing  and  the  Hankel  transform.  We  will  later  use  the  results  derived  in  this  section  to  determine 
the  range  over  which  pressure  field  data  must  be  known  in  order  to  successfully  estimate  the 
plane  wave  reflection  coefficient.  We  will  also  later  use  the  approximate  results  presented  in  this 
section  to  determine  the  effect  of  varying  source-height  during  data  acquisition  on  the  estimate  of 
the  plane  wave  reflection  coefficient.  Section  0.7  studies  the  effect  of  sampling  on  the  Hankel 
transform.  Sampling  issues  arise  both  when  data  to  be  transformed  is  available  only  on  a  discrete 
set  of  points  and  when  the  Hankel  transform  is  computed  numerically.  The  results  from  this  sec¬ 
tion  will  be  used  extensively  in  Chapter  IV. 


f 
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The  addition  of  white  Gaumian  noise  is  often  a  reasonable  model  for  the  accumulated 
effect  of  many  sources  of  corruption  acting  on  a  measured  signal.  Section  n.8  discusses  the 
degradation  introduced  into  the  Hankel  transform  of  a  signal  by  the  addition  of  white  Gaussian 
noise.  It  also  shows  that  sampling  such  a  function  on  a  square  root  grid  can  improve  the  noise 
behavior  of  the  associated  Hankel  transform. 

We  begin  now  by  presenting  common  definitions  of  the  Hankd  transform. 


D.2)  Definition  of  the  Hankel  Transform 

ha  the  literature  a  number  of  different  integral  transforms  ate  referred  to  as  the  Hankel 
transform.1  Three  of  these  are  presented  below: 


1) HTx 

2 ) HT2 


^(r)J  -  f/(r)Jo(pr)rdr  -  F x(p)  Watson  [1966] 

(r)|  *  2irJ/(r)Jo(2'irpr)rdr  “  F2(p)  Bracewell  [1965] 

3)1173  |f(r)|  *  !/(r)J0(f>r)V^dr  m  F3(p)  Bateman  [1953] 

Definitions  (1)  and  (2)  are  only  superficially  different  since  f  2(p)  -  2<irf‘i(2irp).  Definition  (3) 
is  substantially  different  with 


As  we  will  see,  under  definition  (3)  the  Hankel  transform  has  properties  very  similar  to  the 
Fourier  transform.  We  will  use  definition  (1),  never-the-less,  because  of  its  relationship  to  the 
two  dimensional  Fourier  transform. 


1)  S ometimm  them  transform  art  also  referred  to  as  the  ten-order  Hankel  transform.  We  win  not  make 
(hat  disttoaton  in  this  thesis. 


D.3)  The  Hankel  Transform  as  a  Two  Dimensional  Fourier  Transform 


If  we  use  die  definition  of  Wstson 


HT 


H-b 


>Vd(p r)rdr  ■  F(p)  (1) 

then  die  Henkel  transform  is  uqiljr  related  to  the  2  dimensional  Fourier  transform  of  a  circu¬ 
larly  symmetric  function.  [1,2  ]  To  dtow  this  we  write  the  2  dimensional  Fourier  transform  in 
cartenan  coordinates; 


Fc(kxJcf)  m  JLf  ffc(x,y)el(k'x+k*’)dxdy  (2) 

alt  _»_» 

H  fc  (jr  )  is  circularly  symmetric  we  can  unambiguously  define 

fir)  ■/c(*,y)  where  r  *  y/x2+y2  (3) 

Writing  (2)  in  polar  coordinates  we  have: 

»2tr 

MP.4>)  »  ^///(r)e,preo,(e-*Vdrd0  (4) 

ait  o  o 

A  change  of  variables  £  ■*  0— <J»  shows  that: 

-  ^-nfir)e‘^rdrd%  (5) 

aif  oo 

so  Ff  (p,<b)  is  not  a  function  of  <(>.  We  suppress  <J>,  drop  the  subscript,  P ,  and  perform  the  £ 
integration  using 

i  2w 

ofcc)  (6) 

to  see  that  any  radial  slice  of  the  two  dimension  Fourier  transform  of  the  circularly  symmetric 
function  /c  (x  O' )  is  given  by: 

m 

F(p)  -  f/(r)Jo(pr)rdr  (7) 

o 

which  is  the  Hank  el  transform. 


By  considering  the  Hankel  transform  as  the  two  dimensional  Fourier  transform  of  a  circu¬ 
larly  symmetric  function  we  can  also  relate  the  Hankel  transform  to  the  Abel  transform.  The 
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Abel  transform  frequently  arises  in  optics,  seismology  and  other  fields.  In  this  formulation  it  will 
appear  as  the  projection  of  a  two  dimensional  circularly  symmetric  function  onto  its  axis. 

We  begin  by  noting  that  the  slice  of  the  two  dimensional  Fourier  transform  in  polar  form, 
F/»(p,0),  equals  the  slice  of  the  two  dimensional  Fourier  transform  in  cartesian  form,  Fc(p,0), 
since  both  functions  represent  the  same  slice  of  the  two  dimensional  Fourier  transform.  When 
fc(*,y)  t*  circularly  symmetric  then  its  transform  in  polar  form,  (p,<|>),  is  circularly  sym¬ 
metric  and  equal  to  F  (p),  its  Hankel  transform,  as  we  have  shown.  For  this  case  we  can  there¬ 
for  write: 


F(p)  -  Ml  P.0)  =  M  p,0)  =  J  Sfc(x,y)e^dxdy  (8) 

Alt  _X-I 

If  we  perform  the  y  integration  first  we  have 


5  fcix,y)dy 


ei(adx 


(9) 


The  integral  in  y  generates  the  projection  of  fc  (x  ,y  )  onto  the  Jr  axis.  We  define  this  projection 
to  be  p  (x  ).  If  we  use  the  circular  symmetry  of  fc  (x  ,y  )  we  can  rewrite  this  projection  as: 


X 

/>(*)  =  5  f  A*  >y)<*y 
— » 


a»  __ __  * 

f  fc(V*2+y2'0)dy  =2//c(V^V,0)dfy 


Or  in  the  cylindrical  coordinate  system 


(10) 


P(x) 


(11) 


(r)rdr 

Equation  (11)  is  the  Abel  transform  of  f(r).  The  Abel  transform  can  therefor  be  considered  as 


foe  projection  of  a  circularly  symmetric  fc(x,y)  onto  the  x  axis.  Since  the  Hankel  transform 
was  shown  to  be  the  Fourier  transform  of  die  projection  we  see  that 


F( p)  -  (/O’)  j  J  (12) 

This  relationship  was  presented  by  Bracewell.  [3  ]  Implementation  of  the  Hankel  transform 
through  Equation  (12)  is  equivalent  to  the  projection-slice  method  proposed  by  Oppenheim, 


Friak,  and  Martinez.  [4  ] 

Equation  (12)  relates  the  Hankel  transform  and  the  one  dimensional  Fourier  transform 
through  the  Abel  transform.  When  we  consider  only  large  values  of  p,  the  transform  variable, 
an  approximate  relationship  between  the  Hankel  transform  and  the  one  dimensional  Fourier 
transform  can  be  developed  that  does  not  involve  the  Abel  transform.  This  can  be  done  through 
the  asymptotic  form  of  die  Hankel  transform,  which  we  present  in  the  next  section. 

n.4)  The  Asymptotic  Form 

If  the  Hankel  transform  is  not  dominated  for  all  values  of  p  by  the  behavior  of  the  kernel 
near  the  origin  (as  would  be  the  case  for  h(r)/r  for  example  )  then  the  asymptotic  behavior  of 
the  transform  can  be  studied  by  substituting  in  the  asymptotic  form  for  the  Bessel  function.  If 
we  use  the  asymptotic  form  for  the  Bessel  function  presented  by  Lipschitz  [5,6]  1 

•^Vo(x)  =  cos(*-- J)  +  ~sin(a— J-)  -  cos(x--J)-0^-  x>0  (1) 

where  ]8 1  S  1  and  we  keep  only  the  leading  terms  in  x ,  the  Hankel  transform  becomes: 

VTp]> (p)  ~  (r)cos (jpr  )-  —y^rdr  (2) 

If  we  expand  the  cosine  term  Equation  (2)  becomes: 

so  m 

vTp|>  (p)  ~  T7-  //  (,r  )cos  (pr  Y^rdr  4*  sgn  (p)J*/  (r  )sin  (pr  )v7irfr  (3) 

’f  o  0 

The  integrals  in  Equation  (3)  are  the  Fourier  cotine  transform  and  the  Fourier  sine  transform  [8 
].  In  some  cases  this  form  allows  us  to  extend  results  available  for  these  Fourier  transforms  to 
the  Hankel  transform.  When  die  sgn  (p)  term  can  be  ignored,  for  example,  Equation  (3)  sug¬ 
gests  that  asymptotically  ^|p|F(p)  behaves  much  like  the  Fourier  transform  of  ^  |r  |/(r). 
The  sfn(p)  term  on  not  be  ignored  without  further  approximation  when  the  values  of  the  sine 

1)  A  sort  na«  ntaw  is  the  tens  of  an  asymptotic  mrias  with  tha  tarns  Waiting  tana  it  [7 ) 
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and  cosine  transforms  fox  negative  p  effect  the  poative  part  of  the  spectrum.  Such  is  die  case 
when  these  transforms  are  degraded  by  sampling  or  integration  to  a  finite  limit,  for  example.  [9  ] 

Had  we  used  the  definition  of  Bateman  for  the  Hankel  transform,  Equation  (3)  would  have 
appeared  erven  more  like  a  Fourier  transform: 

at  n 

F3(p)  =  -4-  Jf(r)cos(pr)dr  +  J*n(p)j/(r)rin(pr)dr  (4) 

vw  o  0 

Bateman’s  definition  (Equation  H.2.3)  is  more  directly  related  to  the  Fourier  transform  than  die 
definition  of  Watson  (Equation  n.2.2).  Despite  this,  we  use  the  definition  of  Watson  because 
we  with  to  preserve  the  relationship  between  die  Hankel  transform  and  the  2 -dimensional 
Fourier  transform  presented  in  Section  H.3. 
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n.5)  General  Properties  of  the  Hankel  Transform 

A  number  of  properties  for  the  Hankel  transform  are  readily  available  in  the  literature. 
[10, 1, 6, 7]  We  present  some  of  the  more  important  of  these  hare  for  completeness.1 


PROPERTY 

Hr) 

90 

F(p)  -  J/(r)Joipr)rdr 

0 

self  -"inverse 

F(p) 

f(r) 

linearity 

a  F  i(p)+F  2(p) 

scaling 

/(<«•) 

a*  o 

derivative 

vJ/(0 

“P  2f  (p) 

power 

*  90 

//  (r)g*(r)rdr  =  Jf  (p)G*  (p)pd 
0  0 

P 

moment 

with  m.  *  Jr”f(r)dr 

0 

'W  .3,  («!)V  P 

hi  the  remainder  of  this  chapter  we  develop  properties  of  the  Hankel  transform  not  avail¬ 
able  in  the  literature  but  of  considerable  importance  to  the  later  developments  in  this  thesis.  We 
begin  by  determining  the  effect  of  on  die  Hankel  transform  of  a  function  when  it  is  multiplied  by 
a  range  limited  window. 

II.  6)  Windowing  and  the  Hankel  Transform 
a)  An  exact  windowing  expression 

The  definition  of  the  Hankel  transform  has  infinity  as  the  upper  limit  of  integration.  In 
practice  it  is  often  impossible  to  carry  out  the  integration  to  infinity.  This  may  be  because  the 
function  to  be  transformed  is  only  known  out  to  a  finite  range  or  because  the  integration  must  be 

*We  will  ooadder  two  function*  to  bo  equal  if  the  remit  of  ooavoMai  their  difference  with  i  bend-limited 
function  ii  always  am.  This  U  equality  in  the  tame  of  (enartUaed  function*. 


i 

l 
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performed  numerically  and  only  a  finite  number  of  calculations  can  be  made.  Following  the 
convention  used  with  die  Fourier  transform  we  will  write  the  upper  limit  of  integration  as  infin¬ 
ity  but  win  make  the  function  to  be  transformed  zero  beyond  some  finite  upper  limit  by  multi¬ 
plying  by  a  window  of  finite  extent.  [9  ]  In  this  section  we  will  explore  the  degradation  intro¬ 
duced  into  the  Hank  el  transform  of  a  function  by  such  windowing.  The  results  of  this  section 
will  also  find  application  to  the  approximate  evaluation  of  integrals  of  the  form 


9 

SfirVfor)J<ter)rdr  (1) 

0 

which  arise  in  this  thesis  in  connection  with  source-height  effects. 

An  exact  but  cumbersome  expression  for  the  effect  of  windowing  can  be  derived  from  a 
result  presented  by  Bracewell.  [10  ]  If  we  define: 


9 

HP)  m  Jp(r)J(^pr)rdr 
0 
9 

W( p)  *  /w(rVo(pr)rdr 

o 

Then Pw(p),  the Hankel  transform  of  the  product  otp(r)  and  w(r)  is  given  by: 


(2) 


*  »2*  __________ 

PAP)  -  Jp(r)^(r)J^pr)rdr  «  ffp (©  W (V p2+ t2-2p^oos8)^ dd $  (3) 

0  0  0 

We  can  relate  the  Hankel  transform  of  the  windowed  function,  Pw  (p)  to  the  Hankel  transform 
of  the  unwindowed  function,  P(p),  by  carrying  out  the  theta  integration  in  Equation  (3)  to 


obtain: 


»  2* 

with  Hi p,€)  -  i£w (Vp2+42-2p$cose)</e  (4) 

If  H(p,4)  had  the  form  /?( p-£),  then  Equation  (4)  would  be  a  convolution,  reminiscent 
of  the  windowing  result  for  the  Fourier  transform.  [9  ]  By  placing  some  restrictions  of  w  (r)  we 
can  derive  an  approximate  expression  for  the  effect  of  windowing  that  has  the  form  of  a  convo¬ 


lution. 
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b)  Approximation  as  a  convolution 

A  simpler  approximate  expression  describing  the  effect  of  windowing  can  be  derived  by 
using  the  asymptotic  expression  for  the  Hankel  transform: 

■  9 

vTp"|F(p)  *  ^  ff(r)cos(pr)v7dr  +  sgn (p)|/ (r )jin (pr )^7dr 

Gt  (p)  and  Gq(p)  are  a  Fourier  cosine  transform  and  Fourier  sine  transform  respectively.  [8  ] 
The  cosine  transform  and  the  sin  transform  each  have  the  property  that  the  transform  of  a  pro¬ 
duct  is  the  convolution  of  the  transforms.  Using  this,  the  effect  of  windowing  in  die  asymptotic 
formulation  of  Equation  (1): 

*  9 

V^IMp)  ~  —  Sf(r)w(r)cos(pr)V^dr  +  sgn (p) // (r )w (r )sin (pr )^7dr  (2) 

0  0 

can  be  written  as: 

VIpI>*(p)  *  G.(p)MMp)  +  sgn( p)  [G0(p)*WJr(p)]  (3) 

where  we  have  defined: 

3B 

*>(P)  ■  Jw(r)eifrdr  (4) 

—3* 

In  general  Equation  (3)  can  not  be  rewritten  as  the  convolution  of  F(p)  with  a  window 
term  because  of  the  sgn  (p)  term.  However,  if  die  Fourier  transform  of  the  window,  Wp(p),  is 
effectively  confined  to  a  narrow  band  around  p  -  0,  then  for  p  larger  than  this  band  (fi„): 

P  >*w  Gt(p)*WF(p)+sgn{f>)  [o0(p)*W>(p)]  *  G,(p)*W>(p)+  [j*n(p)G0(p)  ]*WV(ptf) 
Combining  Equations  (2),  (3),  and  (5)  we  have: 

''TpIMp)  88  [vTpI*’(p)]*«V(p)  (6) 

which  is  our  asymptotic  result. 

If  IFjr(p)  is  not  negligible  beyond  some  band,  Bm,  then  the  effect  of  windowing  on  the 


-  34  - 


Hankel  transform  can  not  be  put  in  this  simple  form.  For  such  cases  the  exact  result  of  Equation 
(n.6a.4)  must  be  used, 
c)  Resolution  and  leakage 

Given  Equation  (II. 6b. 6)  we  can  address  die  practical  issues  aaodated  with  windowing. 
Am  is  frequently  done  for  the  Fourier  transform,  we  divide  the  issues  associated  with  windowing 
into  two  general  classes.  The  first  we  call  resolution  and  refers  to  die  local  snearing  affected  by 
die  main  lobe  of  the  window.  The  second  we  call  leakage  and  refers  to  the  contribution  of  the 
ride  lobes.  [9, 11  ] 

We  begin  by  expanding  Equation  (I.6b.6)  to  write: 

X  X 

P  >0  ^th  WF( p)«  JW(r)e‘rdr  (1) 

0  -* 

When  p  is  sufficiently  large  (  p  greater  than  some  po)  then  ^  can  be  considered  constant  over 
the  main  lobe  of  WF  (p— £).  For  these  p.  Equation  (1)  can  be  written  approximately  as 

m 

P>P0  ^F*(p)  *  ^pjF(OWF(p-Odi  (2) 

0 

so  that 

Mp)«  J>(0"V(p-0d$  (3) 

0 

Under  this  condition  the  issues  of  resolution  for  the  Hankel  transform  are  the  same  as  those  for 
the  Fourier  transform.  If  we  desire  to  resolve  events  in  the  Hankel  transform  on  the  order  of  5 
then  the  lobe  of  our  window  must  be  less  than  8.  Discuss ons  about  a  variety  of  windows  are 

available  in  the  literature.  [9, 11  ]  For  die  Hanning  class  of  windows,  the  main  lobe  width  is 

where  B  is  die  length  of  the  window.  Our  requirement  for  resolution  of  events  on  the  order  8 
becomes: 

*  >  i  w 

Leakage  is  the  phenomena  we  associate  with  the  ride  lobes.  For  the  purpose  of  this  analysis 
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we  consider  the  lobe  width  to  be  sufficiently  small  that  it  can  be  approximated  by  an  impulse  so 
that  we  ignore  the  smearing  effects  that  we  have  assigned  to  resolution.  We  approximate  Wjr(p) 
as  a  weighted  superposition  of  impulses: 

*V(p)  =  2>/«(p-r,)  (5) 

i 

The  Of  indicate  the  rate  at  which  the  side  lobes  approach  zero.  The  convolution  of  Equation  (2) 
becomes: 

a» 

^(p)  * 

-»  I 

Sa/VF^Xp-r,)  (6) 

i 

When  we  are  concerned  about  the  leakage  due  to  a  singularity,  we  must  consider  die  weighting 
flfVp-T,  which  indicates  the  amount  of  leakage  of  an  event  at  T,  of  strength  1  would  have  at 
p.  Here  the  Hankel  transform  differs  from  the  Fourier  transform  because  of  the  Vp— 7*,-  term 
which  slows  the  decay  rate.  Consequently  for  equivalent  leakage,  the  lobes  of  the  window  must 

fall  by  a  factor  of  —7-  faster  than  that  required  for  equivalent  performance  in  the  Fourier 
VP 

transform.  For  this  reason  we  have  concentrated  on  the  Hanning  window  rather  than  the  Ham¬ 
ming  in  many  of  our  examples. 

By  weighting  the  side  lobe  heights  by  a  factor  of  "^p,  optimal  windows  could  be  designed 
for  Hankel  transforms  in  a  manner  analogous  to  the  Fourier  transform. 


d)  Examples 

In  this  section  we  present  two  examples  of  windowing  and  the  Hankel  transform.  To  each 
we  apply  a  rectangular  window: 

fl  0  <  r  <  4000 

"<r>“\0  4000  <r  W 

which  has  a  length  similar  to  the  range  over  which  data  is  available  in  the  experiment  described 

,flVrJ+2J 

in  Section  (1.6).  The  first  function  we  transform  is  ^  ^  ^  for  which  the  true  Hankel 
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transform  is  given  by  — i  ^  e 1  Figure  n.6d.l  presents  the  log-magnitude  of 

V*2_ p2 

.ftVrJ+(2)J 

— >  0  <  r  <  4000.  As  can  be  seen  in  the  figure,  this  function  decays  almost  four 

Vr2+(2)2 

orders  of  magnitude  over  the  window  length.  Figures  n.6d.2a  and  n.6d.2b  present  the  magni¬ 
tude  and  phase  of  its  computed  transform.  Essentially  no  degradation  due  to  aliasing  is  apparent 

ettVr *+(133)* 

in  the  computed  transform.  Figure  n.6.3  presents  the  magnitude  of  — j  .  .  Over  the 

Vr2+(i33)2 

window  length  this  function  has  decayed  roughly  two  orders  of  magnitude.  Figures  n.6d.4a  and 
IL6d.4b  present  the  magnitude  and  phase  of  its  Hank  el  transform.  The  correct  transform  is 

given  by:  ^  ^  ^el  *2  ^u3)  j^e  magnitude  of  the  correct  transform  should  look  Klee 
glkVri+Q)2 

— j  "■■■  for  0^p£k.  Instead  the  magnitude  of  the  transform  shown  in  Figure  n.6d.4a 

v*2_p3 

shows  considerably  more  degradation  than  that  of  Figure  II.6d.2a.  This  is  due  to  the  fact  that 
this  is  the  transform  of  a  function  which  has  proportionally  more  energy  outside  the  window. 


One  is  tempted  to  assume  the  ripples  apparent  in  Figure  n.6d.4a  are  due  to  leakage  of  the 
j  ^  singularity.  This  is  not  the  source  of  degradation,  however  as  may  be  seen  by  noting 

that  this  same  singularity  is  present  in  the  first  transform  of  Figure  Q.6d.2a  for  which  no  such 
rippling  is  apparent.  The  rippling  is  due  to  the  smearing  of  the  transform  in  Figure  n.6d.3a  over 
its  rapidly  oscillating  phase  term  4»*v*2~p*(133)  wj1jcj1  j,  not  apparent  in  the  magnitude  plot. 
When  the  true  phase  varies  rapidly  over  the  width  of  the  main  lobe  of  the  window  the  effect  of 
smoothing  can  actually  be  to  introduce  rippling  into  the  computed  magnitude. 
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Figure  n.6d.2b  Phase  of  numerically  generated  Hankel  transform  of  function  drown  in  Figure 

n.6d.l 


PHASE 


II.7)  Sampling  and  Aliasing 

It  is  often  necessary  to  approximate  the  integral  in  the  Hankel  transform  by  a  sum.  This 
approximation  may  be  necessary  because  the  function  to  be  transformed  is  known  only  on  a 
discrete  set  of  points  or  because  the  integral  must  be  evaluated  numerically.  The  resulting  sum 
will  be  a  degraded  version  of  the  true  Hankel  transform.  We  will  adopt  the  terminology  of 
Fourier  transforms  and  refer  to  the  replacement  of  the  integral  by  a  sum  as  sampling  and  the 
resulting  degradation  as  aliasing.  In  this  section  we  examine  the  form  that  aliasing  takes  for  the 
Hankel  transform. 


The  discrete  sum  approximation  that  we  will  concentrate  on  is  the  Fourier-Bessel  series. 
We  will  derive  an  expression  that  relates  the  output  of  the  Fourier-Bessel  series  to  the  true 
Hankel  transform.  Because  the  Fourier-Bessel  series  uses  samples  on  a  set  of  points  that  is 
approximately  evenly  spaced,  the  results  we  derive  will  be  approximately  correct  for  any  evenly 
spaced  sampling  scheme. 


We  begin  with  the  formulation  of  the  Fourier-Bessel  series  [7, 6  ]  which  states:1 


. £  , 

0<p<l  F(p)  =  22  J - 777—7 - «fo(*«P) 

« “io  JiKK) 


(i) 


Where  X,  n  =1,2,3,  •  •  •  are  the  ordered  zeros  of  JoC*). 

If  F(p)  =  0  for  p  >  1  then  the  integral  in  the  expression  above  is  just  the  Hankel 
transform  of  F  (p)  evaluated  at  \„ ,  / (X,  )  so  that  the  Hankel  transform,  F  (p),  can  be  expressed 
exactly  as  a  sum: 


»  f  (\  \ 

0<P<1  F(p)  -  2  2  ^TTTTtMKp)  when  F(p)  «  0  for  p  >1  (2) 

■  i  (*») 

When  F  (p)  is  not  truly  bandlimited  to  p  <  1  and/or  the  sum  is  not  carried  out  to  infinity,  Equa¬ 
tion  (2)  is  only  an  approximation  to  the  Hankel  transform.  The  study  of  the  effect  of  finite  N  on 


1)  We  will  call  two  fractions  equal  if  the  Fourier  trasiom  of  hair  difference  has  00  energy  at  say  finite 
frequency.  For  this  reason  wa  need  not  single  out  the  values  of  F  (p)  in  Equation  (1)  at  points  of  disoon* 
Sanity. 
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the  approximation  is  the  study  of  windowing,  covered  in  the  previous  section.  Here  we  consider 
only  the  degradation  that  occurs  because  the  infinite  series  is  used  in  place  of  the  integral. 

Finally,  we  note  that  it  is  because  the  zeros  of  Jq(-x  ),  Kn,  rapidly  approach  nit— ~  that  the 

sampling  above  is  approximately  evenly  spaced. 

To  determine  the  effect  of  approximating  the  Hank  el  transform: 


F(p)  -  ff(r)Jd.pr)rdr 
0 

by  the  Fourier-Bessel  series: 


(3) 


0  <  p  <  1  F( p)  -  2  /(X.V 0(KP) 


(4) 


we  express  F  (p)  in  terms  of  the  correct  transform,  F( p),  by  inverting  (3)  to  write  f(r  )  in  terms 
of  F(p).  We  substitute  this  into  Equation  (4)  to  yield: 


N 


0  <  p  <  1  F(p)  =  2 


'lJUK) 


fFWoiKZKdti 

0 


Interchanging  the  order  of  integration  and  summation  we  have 


^o(^«  P) 


(5) 


/(p)  -  j>(«)r*<p,oerf« 

0 

Where,  following  the  notation  of  Watson  (page  582)  [7  ]  we  define: 


tn(  p,€)-22 
«»1 


Jofr,tyo(K„p) 

jHk) 


(6) 


(7) 


The  study  of  aliasing  for  the  Fourier-Bessel  series  is  the  study  of  T *(p,£).  We  can  obtain 
an  expression  for  7„(p,€)  by  using  an  asymptotic  result  presented  by  Schlafli:1,2  [12  ] 

[  sin  4»(p-{)  sin  4)»(  2-p—{) 


7V(p,5)  - 


dn  tin 


where  Ah  •  (N+—)v 


(8) 


1)  TUs  differs  from  Wstacn’i  presentation  of  Sdilafli'i  result. 

2)  Schlafli  does  not  restrict  the  regioa  of  validity  of  his  result  Watson,  however,  states  that  Sehlaffi's  result 
proceeds  from  e  formula  which  Is  strictly  valid  only  for  0<p+$<2  sad  p*{.  We  will  later  plot  T m(p,0 
to  show  that  the  results  of  this  analyhs  appear  valid  indde  the  region  0<p+£s2  and  approximately  valid 
outside  that  region. 


3 

1 
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As  N  -*  TN(p,i)  approaches  a  weighted  sequence  of  impulses.  We  determine  that  sequence 
here. 

We  begin  our  analysis  of  Equation  (8)  by  first  considering  the  expression: 

tin  Anx 


sin 

..  ,  irx 

Afirx  +  -r- 

4 

tin 


irx 


iin— 


(9) 


Which  equals 


sin  JV  irx  __  irx  .  cos  N irx  irx 
■■■"■  COS  — “  +  ■  —  Sin  -r- 

irx  4  irx  4 


(10) 


sin 


sin 


2  ~  2 

In  Appendix  I  we  show  that  as  AT  -  *  the  first  term  in  (10)  approaches  the  limit: 


2(-l)*8(~  -2*) 

k  * 


(ID 


hi  Appendix  I  we  also  show  that  the  second  term  in  Equation  (10)  approaches  0.  The  limit 
of  Equation  (9)  is  therefore  given  by: 


Ar-»  irx  l  l 

sm  ’ 

2 

Using  Equation  (12)  the  first  term  in  Equation  (8)  can  now  be  seen  to  approach  the  limit: 

sin  AN( p-0 


(12) 


Km 

nn  *(P~£) 


22(-l)*5(p-«-4Jt) 


(13) 


The  second  term  in  Equation  (8)  can  be  put  in  the  form  of  Equation  (9)  by  defining 

y  m  2— p— 4: 

«in  An  (2-p~4)  sin A^y 


2  2 


(14) 


Combining  Equations  (13)  and  (14)  we  have: 

sjnMy(2-p-4) 


&n  "  22(-1)*8(2-p-4-4*) 

^  *(2-p-4)  * 


(15) 
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We  can  determine  7*(p,£)  by  combining  (13)  and  (15):1 

^(p.0  =  ^rjr2(-l)*  [»(P~€-4*)  -  8(2-p-£-4*)]  (16) 

If  our  transform  is  not  severely  aliased  so  that  F  (p)  is  negligible  for  p  >  2  then  substituting 
Equation  (16)  into  Equation  (6)  shows  that: 

0<p<l  /(p)«/FU)^L-[8(p-5)  -8(2-p-S)]$dS  (17) 

which  equals  for  0<p<l  :2 

*  (p)  -  r  (p)  -  (2-p)  (18) 

We  observe  that  the  aliasing  result  most  directly  relates  F (p)  to  V^F(p). 

An  example  of  aliasing  is  presented  in  Figure  Q.7.1  where  we  see  4'v^p  times  the  Hank  el 
transform  of  e-*'2  generated  with  the  Fourier- Bessel  series.  The  figure  displays  die  aliasing 
terms  generated  by  the  impulses  in  Equation  (16).  In  the  region  0<  p<  2  the  figure  mutches 
the  result  indicated  in  Equation  (17)  very  well.  In  the  region  0<  p<  4  the  figure  does  not 
correspond  exactly  to  what  would  be  determined  by  substitution  Equation  (16)  into  Equation  (6) 
indicating  the  limited  validity  of  Schlafli’s  result. 

Figure  H.7.2  shows  a  plot  of  2‘^/p£7'i28(p,£)  0<  p<  10  0<  £<  10.  This  picture  sup¬ 
ports  the  accuracy  of  Equation  (16)  for  7*(p,£)  for  0<  p+£:£  2  and  suggests  that  Equation 
(16)  is  at  least  approximately  correct  over  the  range  of  p  and  £  shown  in  the  figure. 

We  conclude  this  section  with  a  final  example  of  aliasing  for  the  Hankel  transform  that  will 
play  an  important  role  in  the  generation  of  synthetic  data.  Figure  n.7.3  shows  the  function 

/(')  “  "J1  a  Im(r0)  >  0  (19) 

corresponding  to  two  poles,  one  at  r  =  r0  and  the  other  at  r  *  — r0.  This  function  has  the 
known  Hankel  transform:  [8  ] 

1)  Over  the  region  of  validity  for  Schlafli'i  result. 

2)  We  have  included  tha  point  p+£  ”  2,  which  it  not  itrictly  within  tha  interval  qpedfied  by  Waaon. 


MAGNITUDE 


t 


Figure  11.7. 1  A'^pHT  je  _SrI  |  generated  numerically 
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f(p)  =  ~Hil)M 


(20) 


Vj  —  lr 

Asymptotically  ff^^oP)55  , )  e  *  e  80  that  its  magnitude  should  appear 


"irr0p 

V2  1 

■>.  Consequently  the  magnitude  of  F( p)  should  appear  smooth  and  decay  as  —r*.  hi 
/irrop  vn 


P 

.  Rapid 


Figure  II. 7. 4  we  see  the  numerically  computed  transform  using  the  samples  / 


K 


2048 


oscillations  are  apparent  which  are  not  present  in  the  magnitude  of  the  correct  transform.  The 
source  of  these  oscillations  is  aliasing,  which  can  be  seen  by  using  Equation  (18)  to  approximate 
the  numerically  computed  transform: 


F(p)  s  ^5  _  v4096-p  r~r~ >„(*096-p) 

V-nr0p  Virr0(4096-p) 


(21) 


which  equals 


^(p)  ~  T7- 

vp  Virr0 


1 


eirop  -Ae~l,ofi 


(22) 


with  A  3  This  can  be  rewritten 


^(P)  88  J'r *  *  [(1  -A)e‘roP  -r  2iA  sin(r0p)] 


(23) 


p  Vwr0 

where  the  beating  caused  by  the  aliasing  is  apparent  in  the  sin(rop)  term.  The  aliased  output 

displays  the  form  of  the  —7—  decay  term  times  an  extremely  degraded  estimate  of  v'pF(p). 

Vp 

When  the  transform  decays  only  at  a  rate  of  of  -y-,  this  example  shows  that  severe  degradation 
due  to  aliasing  can  be  expected. 


11.8)  The  Effect  of  Additive  White  Gaussian  Noise  on  the  Hankel  Transform 
a)  Statement  of  the  Effect  of  Additive  White  Gaussian  Noise  on  the  Hankel  Transform 

ha  practice  it  is  seldom  posable  to  know  exactly  the  function  whose  transform  we  desire. 
Frequently  it  is  ponible  to  model  the  uncertainty  about  the  input  function  by  assuming  that  the 
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Figure  II.7.4  The  log  magnitude  of  the  numerically  computed  Hankel  transform  of  the  function 
shown  in  Figure  n.7.3 
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errors  associated  with  each  sample  are  random  and  uncorrelated  from  point  to  point.  Since  the 
combined  effect  of  many  random  factors  can  often  be  modeled  with  a  Gaussian  distribution  (by 
invoking  the  central  limit  theorem  [13  ])  the  assumption  is  often  added  that  the  distribution  of 
error  around  each  point  is  Gaussian.  This  model  of  the  uncertainty  corresponds  to  additive  white 
(the  uncorrelated  assumption)  Gaussian  noise.  We  assume  the  mean  of  the  noise  process  is  zero 
so  that  the  expectation  of  the  noisy  input  signal  is  the  true  input.  If  we  further  assume  that  the 
variance  of  the  noise  process  is  not  a  function  of  the  input  sample  number  then  this  Gaussian 
noise  process  is  stationary. 

In  this  section  we  explore  the  effect  of  such  uncertainty  on  the  Hankel  transform  of  the 
input  function.  Since  the  Hankel  transform  is  a  linear  operator  and  the  noise  process  has  zero 
mean,  the  effect  of  the  noise  will  be  to  introduce  a  variance  in  the  output  of  the  Hankel 
transform  proportional  to  the  noise  power  but  the  expected  output  will  not  be  corrupted.  [14  ]  In 
this  section  we  first  show  that  unlike  the  Fourier  transform  the  variance  of  the  Hankel  transform 
of  stationary  white  Gaussian  noise  is  not  stationary,  but  instead  concentrates  power  near  the  ori¬ 
gin.  This  result  is  important  because  frequently  the  Hankel  transform  is  used  in  place  of  the  two 
dimensional  Fourier  transform  in  problems  with  an  underlying  circular  symmetry.  Because  of 
this  property,  a  slice  of  the  the  two  dimensional  Fourier  transform  of  noisy  measurements  made 
over  a  two  dimensional  grid  of  a  circularly  symmetric  field  will  differ  from  the  Hankel  transform 
of  a  slice  of  that  field. 

In  Section  (b)  we  will  show  that  if  /  ( VT )  is  a  stationary  white  Gaussian  noise  process  then 
F  (v  p)  will  also  be  stationary  white  Gaussian  noise.  This  result  implies  that  if  the  input  function 
is  sampled  on  a  v7  grid  and  each  sample  is  independently  corrupted  by  (zero  mean)  Gaussian 
noise  that  does  not  depend  on  the  sample  number,  then  samples  of  the  Hankel  transform  on  a 
vri>  grid  will  be  indepencL.  Jy  corrupted  by  Gaussian  noise  and  the  amount  of  corruption  will 
not  depend  on  the  value  of  p.  On  these  grids  each  sample  represents  the  same  area  of  the  under¬ 
lying  two  dimensional  circularly  symmetric  function  and  the  noise  properties  of  the  Hankel 
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transform  are  equivalent  to  the  noise  properties  of  the  underlying  two  dimensional  Fourier 
transform. 


To  show  that  the  Hankel  transform  concentrates  noise  power  near  the  origin  we  first  write: 

(p)  ■  /  [/O')-*-" (»-)]<Jo(pO«<r  (1) 

Where  we  have  introduced  the  limit  of  integration,  B ,  to  insure  convergence.  n(r)  is  stationary 
white  Gaussian  noise  with  variance  N0.  The  variance  of  Fa  is  given  by 


Mr 


[f.(p)|  - 

|/,(p)-£(fr.(p))l*] 
l/»  (r  Vo(pr  Wr  P  j 

(<*)"  (3)  \j  o(Po  V  0(pp)apd  ad  p 

a  a 

ffN<fi(a -PVo(pa  Vo(pP)a&dad  p 


£ 

f  I 


0  0 


=  iVo^(pa)oJ</a 

For  p  =  0  Equation  (2)  above  shows  that 


VA/?[f,(0)]  -  /Vo/ujda  = 


N  oB3 
3 


When  p*0 


VAX  (p)  1  =  Ar0/7^(pa)aJdo  = 

1  J  0  P  0 


In  units  of  normalized  frequency  v  ■  p/B 


(2) 


(3) 

(4) 


KAi?  [f,  (p)  1  -  N0fj$ (po)alda  =  (Mi  (5) 

1  J  0  Vo  o 

which  is  plotted  in  Figure  H.Sa.l.  As  can  be  seen  this  function  decays  rapidly  with  v  so  that  the 
Hankel  transform  concentrates  noise  power  near  the  origin.  We  can  explain  why  this  noise  pro¬ 
perty  of  die  Hankel  transform  differs  from  that  for  the  Fourier  transform  by  ennaideriBg  the 
underlying  two  dimensional  circularly  symmetric  Fourier  transform  represented  by  the  Hankel 
transform. 


We  recall  from  Section  (11.3)  that  the  Hankel  transform  of  a  function,  /  (r ),  corresponds  to 
a  slice  of  the  two  dimensional  Fourier  transform  of  the  function  / (r  ,8)  made  by  sweeping  f  (r) 
around  the  origin  in  two  dimensions  (so  that  /(r,  0)  =  /(r)  for  all  (t).  When  we  generate  the 
Hankel  transform  of  the  noisy  input,  /(r)  +  n(r),  we  obtain  a  slice  of  the  two  dimensional 
Fourier  transform  of  /(r)+«(r)  swept  around  the  origin.  The  result  is  very  different  from 
sweeping  /(r)  around  the  origin  and  then  adding  SWGN  (stationary  white  Gaussian  noise)  in 
two  dimensions.  In  the  first  case  the  noise  field  is  circularly  symmetric,  in  the  second  case  it  is 
not.  It  is  the  symmetry  in  the  underlying  noise  field  implied  by  the  Hankel  transform  that  causes 
the  concentration  of  noise  power  near  the  origin. 

We  will  now  show  that  this  behavior  of  the  Hankel  transform  with  respect  to  noise  can  be 
averted  by  changing  to  a  v7  coordinate  system  for  the  input  and  a  Vp  coordinate  system  for  the 
output.  Samples  evenly  spaced  in  these  square  root  coordinate  systems  have  the  property  that  the 
distance  between  any  two  samples  always  represents  the  same  area  of  the  underlying  two  dimen¬ 
sional  (circularly  symmetric)  function.  Each  noisy  sample  of  the  function  and  its  Hankel 
transform  represents  the  same  amount  of  area  in  the  underlying  two  dimensional  spaces.  Conse¬ 
quently  the  noise  properties  are  equivalent  to  those  associated  with  samples  evenly  space  on  a 
cartesian  grid  (associated  with  the  two  dimensional  Fourier  transform). 

b)  Proof  that  if  f  ( Vr  )  is  stationary  white  Gaussian  noise  then  F  (^j>)  will  also  be  stationary  white 
Gaussian  noise,  where  F( p)  is  the  Hankel  transform  of  f(r) 

The  proof  that  if  /  (V'r  )  is  stationary  white  Gaussian  noise  (SWGN)  then  F  (Vp)  is  also 

SWGN,  consists  of  three  parts  and  a  conclusion.  First  we  will  show  that  for  the  integral 

* 

transform  defined  as  F2(p)  ■  ff(r)Jo( pr)vprdr  that  /(r)  is  SWGN  if  and  only  if  F2(p)  is 

o 

SWGN.  Second  we  will  use  this  to  show  that  ^rf  (r )  is  SWGN  if  and  only  if  v'p F (p)  is  SWGN. 
Finally  we  show  that  if  v7 f(r)  is  SWGN  then  /  (Vr  )  must  be  SWGN. 

i)  Proof  that  Fj(p)  is  stationary  white  Gaussian  noise  if  and  only  if  /(r)  is  stationary  white  Gaus- 
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sian  noise,  where  F2(p)  the  Hankel  transform  of  f(r)  as  defined  by  Bateman. 

We  begin  by  recalling  that  the  Hankel  transform  defined  by  Bateman  [8  ]  and  presented  in 

x 

Section  (II.2)  is  given  by  F2(p)  ■  Jf(r)j0(pr)^prdr.  Because  this  operator  is  linear,  F2(p)  is 

o 

Gaussian.  For  the  same  reason  if  the  mean  of  /(r)  =  0  then  the  mean  of  F2(p)  is  also  0.  We 
need  only  show,  then,  that 


Now 


E 


^(PiJFrfpj)]  =  CS(p,-p2) 

[fi(pi)F:(pi)  J  » 
f*  * 

JV  ('  )>o(Pl r  )V^drJf  (f )J P(pt0  £ 

0  0 

9  » 

/JjVo8(i- -?Vo(p1r)y0(p2?)Vp,p2r^rd{ 

0  0 

m 

N  o^PiPs  Jj  o(p»CVo(P2f  )%d  £ 


(1) 


(2) 


Pi 

=  No8(p1-p2) 

This  proves  that  /(r)  SWGN  implies  that  F2(p)  is  SWGN.  The  converse  is  proven  by  noting 
that  the  above  integral  transform  is  its  own  inverse  so  that  the  same  argument  applies,  replacing 
/(r)  with  F 2(p). 

ii)  Proof  that  ^pF  (p)  is  stationary  white  Gaussian  noise  if  and  only  if  '^rf  (r)  is  stationary  white 
Gaussian  noise,  where  F(p)  is  the  Hankel  transform  of  /(r) 

From  (i)  we  know  that  if  Vr/  (r )  is  SWGN  then 


J (r )P o(pr y^prdr  =  'Spff(r)Jfor)rdr  -  V^F(p) 


(3) 


is  SWGN.  Since  the  Hankel  transform,  Jf  (r)J o(pr)«fr  is  its  own  inverse  the  converse  follows  by 

0 

starting  with  v9/(r). 
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ill)  Proof  that  'Srf  (r)  is  stationary  white  Gaussian  noise  if  and  only  iffC^r)  it  stationary  white 
Gaussian  noise 

&(rl“r2> 


First  we  show  that  if  E  j/  {j\)f(r$  j  =  S(ri-r2)  then  E  [/(/M)/(/[r2])  j  =  —j 
This  assumes  that  /(r)  and  f-1(r)  exist  in  the  region  of  interest. 


Proof: 


fOlr  1)  = 

s 


Using  Equation  (4)  we  can  write  E  j/  {l[r\])f  (/ [r2])  J  as: 

[/(/[ri])/(/[r2])]  - 


*1 
W  » 


J*  /8(a1-tt1) 


B[r1-f~1(a1))8{r2-/~>(tt2))4ai4a;i 
8[ri~ t~*(ai)]8[r:»— /~*(a;)]4ai4a; 


f  a[r«~rlfr)l6trs-rl(a)] 
-»  /’[/'‘(a)] 


If  we  define  £  ■  /_1(«)  «i  that  a  =  l  (i)  and  da  =  l (£)d £  then 


*(rx-rf) 

i(ri) 


We  now  show  that/O^ )  SWGN  iff  E  j/(»'i)/(r2)  J 


no 


r\ 


(4) 


(5) 


(6) 


We  apply  the  result  of  the  first  part  of  this  section  to  our  special  case  by  defining  l(r)  ■  r2 
so  that  i(r)=*2rdr  With  this  definition  we  see  that  if  some  p(r)  is  SWGN  then  p(r2)  will  have  a 


variance  proportional  to  ~ .  This  implies  that  if  /  (V7 )  is  SWGN  then  /(r)  has  a  vari- 
*(ri~r2>  _ .v..  «  vl  &(r»~ra) 


once  proportional  to  — ^ — —  .  Finally  we  note  that  if  E  |/‘(r1)/(ra)  j  -  —  ^  —  then 
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£[V^(ri)V^(rj)]  *  8(rj-r2)  therefor  If  / (v7 )  is  SWGN  then  v7/(r)  is  SWGN.  The 
converse  follows  from  the  same  steps  in  reverse. 

iv)  Conclusion 

In  (i)  and  (ii)  we  showed  that  if  Vr/ (r)  is  SWGN  then  V/p F (p)  will  also  be  SWGN.  In 

(iii)  we  showed  that  'Srf  ( r )  is  SWGN  if  and  only  if  fC^r  )  is  SWGN.  The  argument  of  (iii)  and 

(iv)  also  showed  that  VpF(p)  is  SWGN  if  and  only  if  F(V p)  is  SWGN.  Thus  we  have  shown 
that  if  / (v7 )  is  SWGN  then  F (Vp)  must  also  be  SWGN. 


II.9)  Summary 

hi  this  chapter  we  have  developed  the  properties  of  the  Hank  el  transform  in  general.  By 
themselves  these  properties  might  be  simply  interesting  but  together  with  the  ability  to  numeri¬ 
cally  perform  the  Hankel  transform,  they  become  tools  for  scientific  research.  In  the  next  chapter 
we  will  discuss  methods  for  numerically  performing  the  Hankel  transform. 
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Chapter  III: 

Computing  the  Hankel  Transform 

m.l)  Overview 

Because  of  the  importance  of  the  Hankel  transform  there  presently  exist  many  algorithms 
for  its  evaluation.  As  yet,  however,  there  is  no  generally  accepted  algorithm  analogous  to  the 
FFT  for  the  Fourier  transform.  This  chapter  presents  a  survey  of  numerical  Hankel  transform 
tedmiques.  It  does  not  indude  all  the  published  algorithms  but  does  represent  the  major  classes. 
We  begin  in  section  HI. 2  by  discussing  perhaps  the  oldest  and  best  understood  algorithm,  the 
Fourier-Bessd  series,  that  we  used  to  derive  our  aliasing  result  of  Chapter  n.  Section  m.3 
presents  the  backsmear  method  and  an  example  of  an  efficient  class  of  realizations.  Section  m.  4 
discusses  the  asymptotic  transform  method  as  it  has  been  proposed  in  the  literature  and  presents 
new  results  that  can  be  used  to  improve  this  method  or  to  asses  the  error  in  the  standard  realiza¬ 
tion  from  the  output  transform  alone.  Section  TII.S  discusses  a  common  combined  transform 
method  and  presents  a  caution  about  its  use.  Section  IB.  6  presents  a  convolutional  method  fre¬ 
quently  used  for  electromagnetic  problems  which  require  the  transform  of  smooth  functions  of 
limited  extent.  Section  III. 7  discusses  the  projection-slice  method.  In  that  section  we  develop  a 
fast  algorithm  for  generating  the  projections  (shown  to  be  an  Abel  transform),  which  itself  has 
wide  applications.  [1 }  When  it  is  followed  by  an  FFT  the  result  is  an  efficient  Hankel  transform. 
Both  the  Abel  and  Hankel  transform  algorithms  are  illustrated  with  examples. 


m.2)  The  Favrier-Beuel  Series 


Probably  the  first  proposed  method  for  evaluating  the  Hankel  transform  is  the  Fourier- 
Bessel  series,  which  was  discussed  in  the  aliasing  section.  The  Fourier-Bessel  series  is  summar¬ 
ized  by  the  identity:  [2, 3  ] 


0<P<1  f(p)-2i 
»-io 


(i) 
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Where  X„  n  =1,2,3,  *  -  *  are  the  ordered  zeros  of  Jq(x). 

This  relation  can  be  used  to  obtain  the  Hankel  transform  of  a  band  limited  function,  /(r), 
defined  by 


F  (p)  *  //  (r  )J o(Pr  )rdr  (2) 

o 

where  F(p)  =  0  for  p>l  by  noting  that  when  /'(p)=0  p>l: 

*  1 

/(X„)  =  /F(pyo(pX„)p^P  =  J>(p)J0(pX„)p<fp  (3) 

0  0 

Substituting  into  Equation  (1)  we  have: 

0<p<l  F(  p)  =  2  2  757~/0(^  P)  (*) 

which  is  in  the  form  of  a  sum  as  we  desired.  If  F  (p)  is  bandiimited  to  p <B ,  instead  of  p<l. 
Equation  (4)  can  be  modified  by  a  change  of  variables.  The  resulting  general  form  states  that 
when  F  (p)  =  0  for  p>fl  then 


o<P<a  f(p)  *  i  i (5) 

The  properties  of  the  Fourier-Bessel  series  have  been  extensively  studied.  [2, 3  ]  As  a  numerical 
algorithm  for  implementing  the  Hankel  transform  it  is  usually  appropriate  only  when  a  few 
values  of  F  (p)  are  required  or  when  computation  cost  is  less  important  than  careful  control  over 
the  errors.  This  is  because  the  Fourier-Bessel  series  requires  than  a  new  sum  be  computed  for 
every  value  of  p.  Further,  it  requires  that  the  function  to  be  transformed  be  available  on  the 


nonuniform  grid, 


Another  algorithm  that  is  appropriate  when  only  a  few  values  of  F  (p)  are  required,  but 
which  accepts  input  values  on  an  even  grid,  is  the  backsmear  method.  We  present  the  backsmear 
method  in  the  next  section. 
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ni.3)  The  Backsmear  Method 

The  backsmear  method  exploits  the  efficiency  of  the  FFT  to  provide  an  efficient  algorithm 
when  only  a  few  output  values  of  the  transform  are  desired.  This  algorithm  can  be  derived  by 
replacing  the  Bessel  function  with  its  integral  representation  the  Hankel  transform  can  be  written 
as  the  two  dimensional  integral: 

*  lv 

F(p)  *  ff(r)JQ(pr)rdr  = 

o  m  o 

30 

U  we  define  G( x)  m  J*/(r )rein dr  and  note  that  oos6  is  even,  the  transform  becomes: 

o 

F(p)  38  —/G (pcos0)</0  =  J*  |G(pcos8)-rG(-pcos8)  jd0  (2) 

G  (x  )  can  be  evaluated  efficiently  on  an  even  grid  using  an  FFT.  [4  ]  The  integration  called 
for  in  Equation  (2),  however,  must  be  performed  for  every  desired  value  of  p.  Further,  since  the 
integration  of  Equation  (2)  is  in  0,  some  interpolation  scheme  must  be  used  to  generate 
G  (pcos0)  on  the  quadrature  points  required  by  the  numerical  integrator. 

The  backsmear  method  is  appropriate  when  the  transform  is  required  at  only  a  few  values 
of  p.  When  this  is  the  case  and  when  G  (x)  is  slowly  varying  so  that  the  numerical  integration 
which  must  be  carried  out  upon  it  can  be  done  efficiently,  then  the  backsmear  method  is  effi¬ 
cient.  This  is  particularly  true  if  the  interpolation  scheme  used  permits  a  quadrature  formula  to 
be  developed  for  the  numerical  integration. 

One  such  development  is  presented  below  for  the  case  of  linear  interpolation  between  the 
potato  of  L(x)  ■  G(x)  +  G(-Jt). 

We  seek 

jr 

1  * 

F( p)  -  ~/L(pco«8)rfe  (3) 

It  0 

and  have  available  only  L  {~~)  from  the  output  of  the  FFT.  We  approximate  L(x)  by  linearly 

NT 


ff(r)e,prtm*rdr 


d  0 


(1) 
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interpolating  L  so  that 


f- 


0  =  £»-  2  £.,(*) 
il  Jm  0 


where 


I(J_)  +  (x — i-)(AT) 
y  NT  v  NTJ(  ' 

L  ^  ) 

'  NT  '  NTJ 

0  otherwise 

V*) 

We  can  integrate  L  (pcosB)  as  follows 


« 

2 


Tf- 


T-f 


j+l 

■<X  <*■  — 

AT  AT 


J"l(pcos8)d0  =  /  2  iy(pcos0)</0  =  2  J*£/(pcos0)d0 
o  o  ]•  l  /-i  o 

The  Lj  terms  can  be  integrated  term  by  term: 

A  ^-1  i 

■  .Q  hw>  -  4(  $■’  •  t(w>}  * 


t*r 


where 


(4) 


(5) 


(6) 


(7) 


c6s-lx 


(cos-1. 

P 


x  |x  |Sl 
X  121 


Evaluating  this  integral  and  substituting  into  Equation  (6)  we  have 


(8) 


F  (—)  = 
KNT ’ 


2 

J-  0 


cos-1(-^-)  -  cos-1(J 


^L(  AT  )  ~L^NT^ 


II 


(9) 


+  |[vpzj3 .  vpq-^s]  [t(JS,  -  njL) 

When  k  is  mall  or  P(~r)  is  desired  for  only  a  few  values  of  k ,  this  is  an  efficient  peo 


eedure.  When  k  is  large  and  (~~r)  is  required  for  many  values  this  is  a  slow  method,  how* 

NT 


ever.  In  the  next  section  we  present  a  fast  algorithm  for  the  evaluation  of  the  Hank  el  transform 
at  large  values  of  the  transform  variable,  p. 


in. 4;  The  Asymptotic  transform 


In  the  seismic  community  the  Hankel  transfoim  is  commonly  evaluated  approximately  by 
replacing  the  Bessel  function  with  its  asymptotic  form  so  that  the  resulting  integral  looks  like  a 
Fourier  transform.  [5,6  ]  The  result  will  usually  asymptotically  approach  the  Hankel  transform, 

though  pathological  functions  can  be  found  for  which  this  is  not  the  case  (  for  example). 

The  main  disadvantages  of  this  method  are  that  it  is  almost  never  suitable  for  small  values  of  p 
and  that  the  error  induced  by  the  substitution  is  determined  in  part  by  the  function  being 
transformed.  The  main  advantage  is  that  the  resulting  integral  looks  very  much  like  a  Fourier 
transform  and  can  be  evaluated  efficiently  with  an  FFT. 

An  expression  for  the  Bessel  function,  suitable  for  an  asymptotic  expansion  is  provided  by 
Lipschitz:1  [7  ] 

for  x  >0  Jq(x)  =  -^rL  Jcoifi— ~)  +  2-)  -  |-)  +  j  where  I®I<1  0) 

In  the  literature  describing  this  technique  [5  ]  only  the  leading  term  is  retained,  providing: 

F(p)  ==  //  (r  )cos(p r  —  ^-)'vrrdr  when  p  >  0  (2) 

vpno  4 

This  is  written  as  a  Fourier  transform  by  noting  that 

cos(pr-~)  =  ^(coi(pr)  +  si»t(pr))  (3) 

hence 


J/(r  )Vr  cos(pr)dr  +  J/(r)Vr  sin(pr)dr 
o  0 


when  p  >  0 


(4) 


This  expression  can  be  evaluated  with  an  FFT. 

We  can  write  the  error  associated  with  the  transform  technique  as 


«(p)  ■  //('■)»’ 


Mw)  “  yfi-'  cos(pr-y) 
virpr  # 


jff(r)N(pr)dr 


(5) 


1)  A  mart  reamt  expansion  with  the  same  leading  behavior  can  be  found  in  Wilson.  [2  ) 


i 

't 
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with  H  (pr)  »  pr  L’c(pr)  —  r?  cos(pr  —  ~)  ] .  //(.*)  is  plotted  La  Figure  III. 3.1.  It  looks 
L  v  1TPr  4  j 

very  much  like  a  slowly  decaying  sinusoid.  e(p)  decays  with  p  both  because  tf  (pr)  decays 
(slowly)  with  p  and  because  of  the  of  the  1/p  term.  Because  of  the  sinusoidal  behavior  of  H (pr ) 
the  error  term  of  Equation  (5)  can  be  large  at  those  frequencies  where /(r)  has  a  lot  of  energy. 

By  using  the  expression  of  Lipschitz  we  can  develop  an  asymptotic  Kankel  transform  with 
an  error  term  that  decays  much  faster  than  the  error  term  for  the  conventional  asymptotic 
method.  This  expression  also  makes  it  relatively  easy  to  judge  the  validity  of  the  conventional 
asymptotic  transform  after  it  has  been  performed  because  it  presents  the  leading  error  terms  of 
that  transform  as  a  function  of  the  transform  itself. 

We  substitute  the  expression  of  (1)  for  7o(pO  in  the  Hankei  transform  to  see  that: 


F(p)  = 


V2 


128 


J/(r)v7coj(pr-j)aV  -r  j-f^P-s:/i(pr-^-)dr 

hf^yz~cos  “  2t^dr 

r  n  4 


8p  o 


(6) 


If  we  use 


then  Equation  (6)  becomes 


cos(pr-~)  =  [cos(pr)  +  sin(pr)  j 
sin(pr- j)  =  [sin(pr)  -  cos(pr)] 


If  we  now  define 


(V 


(8) 


Mp)  “  J/(r)v7c°s(pr)dr 

0 

* 

F,(p)  ■  J/(r)v7sin(pr)</r 
0 


(9) 
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and  note  that 


then 


T-/ ^P~cos  (PrVr  =  ~ff(r)V7sut(t>r)dr  =  -F,  (p) 

n  v  T  a 


*  //  \  * 

-rJv^sin^dr  =  J/(r)v™M(pr)dr  =  Fc(p) 

"P  0  r  0 

=  -F.(P) 

dp  0  ^ 

77.f,^tt'",,(Pr)</r  =  -F,(p) 

®P  0  r 


(10) 


1 


r (p)  *  ^7  [mp)+f,(p)  +  ^-J>,(?)d?  +  +  jg^JJp.ttddW* 

P  t 


+  -i 

128p-  -  ■  Mp'j 

We  can  combine  terms  above  by  defining 


(ii) 


Fi(p)  «  Fe(p)  +  F,(p) 


so  Equation  (11)  becomes 

...  1 


t 


F (p)  -  ^7[Fi(p)  +  ittW  +  T^JJ>I«1)d&dc  + 


(12) 


(13) 


The  indefinite  integrals  above  can  be  converted  to  definite  integrals  by  integrating  from  0  to 
infinity  and  adding  on  the  values  necessary  to  force  the  resulting  equation  to  match  Equation  (8) 
at  p  =  0: 


p  p  c 

pMey*  -  /FiftW  +  Cl  with  C,  -  J^dr 

a  n  vr 


f  (  p 

JJ>i(fa)dMt  =  / 

o 


ffiCWi  +  C,  d5  +  C2  With  C2  -  j^dr 
.0  o'" 


(14) 


(15) 


Performing  the  integrations  over  the  constants  and  replacing  the  iterated  integral  we  are  left  with 
the  expression 


F(p)  *  ^r- 

v*p 


F«(P)  +  f 


4\|>  i(f)4« 


25 

Ta 


+  J-i 

128  p 


/0»-©Fi(f)4«  4-  C2 


64p3o  r** 


(16) 


The  first  term  can  be  recognized  as  the  usual  asymptotic  expression  for  the  Hankel  transform. 
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The  second  and  third  terms  are  corrections  to  this  expression  which  can  in  principle  be  deter¬ 
mined  from  it.  The  final  term  is  the  remaining  error  term  but  which  is  considerably  smaller  than 
the  term  associated  with  the  usual  asymptotic  expression. 

If  it  were  desired  to  calculate  the  second  and  third  terms  directly  as  a  Fourier  transform, 
the  same  procedure  could  be  used  as  will  be  described  in  Section  (QI.7)  where  it  is  applied  to  cal¬ 
culating  the  Abel  transform. 

Equation  (16)  can  also  be  interpreted  as  providing  first  and  second  order  error  estimates 
for  the  classical  asymptotic  transform.  These  estimates  allow  the  error  associated  with  the  classi¬ 
cal  method  to  be  estimated  (to  first  and  second  order)  from  a  knowledge  of  only  the  resulting 
transform.  Such  estimates  permit  one  to  interpret  the  result  of  an  asymptotic  transform  with  a 
questionable  appearance. 

I1I.5)  A  Combined  Transform  Method 

The  Fourier-Bessel  series  and  the  backsmear  methods  both  permit  the  calculation  of  the 
Hankel  transform  on  any  output  grid.  The  computational  cost  of  each  of  these  methods  increases 
linearly  with  the  number  of  points  computed.  The  asymptotic  method  is  fast  and  usually  gen¬ 
erates  good  estimates  of  the  Hankel  transform  when  p  is  large.  A  combined  scheme  is  possible 
which  uses  a  slow  method  such  as  the  Fourier-Bessel  series  or  the  backsmear  to  compute  the 
Hankel  transform  point  by  point  for  low  values  of  p  and  which  switches  to  the  asymptotic 
transform  for  large  values  of  p.  Such  a  method  has  been  proposed  in  the  literature.  [8  ] 

The  main  issue  with  such  a  scheme  is  the  point  at  which  the  algorithm  should  cease  comput¬ 
ing  the  transform  point  by  point  with  the  slow  method  and  begin  to  accept  the  output  of  the 
asymptotic  transform.  At  present  there  is  no  reliable  method  for  doing  this.  It  has  been  sug¬ 
gested  in  the  literature  that  the  transition  be  made  at  the  first  point  for  which  the  slow  algorithm 
produces  a  transform  value  which  matches  the  value  of  the  asymptotic  transform  within  a  speci¬ 
fied  tolerance.  [8  }  This  scheme  for  switching  to  the  asymptotic  transform  would  be  entirely  ade- 
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quate  if  the  error  made  by  the  asymptotic  algorithm  were  monotonically  decreasing  in  the 
transform  variable,  p.  This  error  is  not  always  monotonically  decreasing,  however,  as  we  will 
illustrate  by  constructing  a  function  for  which  it  is  not.  The  existence  of  such  a  function  indicates 
that  the  switching  method  proposed  in  the  literature  [8  ]  will  not  always  work.  In  fact,  for  the 
function  we  construct,  the  error  made  by  the  asymptotic  method  will  be  zero  at  a  point  we 
specify.  The  switching  scheme  proposed  would  begin  to  accept  the  asymptotic  transform  before 
this  point  even  though  the  erroT  made  by  the  asymptotic  transform  beyond  this  point  might  be 
greater  than  the  specified  tolerance. 

As  in  Section  (m.4)  we  write  the  error  associated  with  the  asymptotic  transform  as: 


«(p)  =  —Jf(r)H(pr)dr 
P  0 


(1) 


where  H  (pr)  *  pr 


[/oCpO  “  -cos  (pr— y) 

virpr  4 


and  was  plotted  in  Figure  HI.3.1.  lit 


order  to  construct  a  function,  /(r),  such  that  e(p)  is  not  monotonically  decreasing  we  first 
choose  some  small  8  and  set 


/(') 


1  0<  r<  8 

s 

— ~SH  (Por)dr 
Po  o 


-Mff(Por)dr 

Po  B 


8<  r<  « 


For  this  /(r)  the  error  made  by  the  asymptotic  transform  is: 


(2) 


m 

«(p)  “  ^SHi»r  )f(r)dr 

(3) 

b  7-/#  (poO^  « 

-  i/»(pr)dr-- - ~fn(pr)dr 

P  0  1  r  .  .  P  8 

Po  8 

When  p  «  po  the  error,  c(po).  equals  zero.  In  general  for  p  >  po  it  will  not.  For  this  con¬ 
structed  function  the  error  is  not  monotonically  decreasing  and  the  switching  procedure  described 


<  ) 
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in  the  literature  will  not  work. 

Until  an  adequate  switching  method  has  been  found  or  until  the  class  of  functions  for  which 
the  proposed  switching  algorithm  will  work  has  been  well  defined,  there  can  be  no  guarantee 
that  the  combined  algorithm  will  work  properly  and  this  method  should  be  used  with  caution. 

HI.  6)  A  Convolutional  Method 

In  this  section  we  describe  a  method  for  computing  the  Hankel  transform  that  puts  the 
Hankel  transform  in  the  form  of  a  convolution  by  transforming  to  an  exponential  grid.  When 
this  grid  does  not  involve  an  extraordinary  number  of  points  to  adequately  represent  the  func¬ 
tion,  the  Hankel  transform  can  be  efficiently  evaluated  with  an  FFT.  We  will  discuss  the  presen¬ 
tation  of  this  method  by  Siegman  [9  ],  though  other  presentations  are  available  in  the  literature. 
[10] 

Siegman  converts  the  Hankel  transform  into  a  convolution  by  sampling  the  function  on  an 
exponential  grid.  He  begins  with  the  Hankel  transform  integral 

X 

F(p)-J/(rVo(pr)rdr  (1) 

0 

After  the  following  definitions: 

P  “  Po*"  *“*»(■£■) 

Po 

r  m  i  ~  fn(~-) 
ro 

Equation  (1)  becomes: 

F(PoO  - 

a 

-  f  f  \ro  *  hV o(Po',o* ,+t )rfc  ^  $ 

— » 

-  r$f  [tf2V(r0«5)]/o(lV(K<,’+<VC 

X 

-  'of  [e_to/('o*"«)]jo(Po'o*(r~8))<fa 


(3) 
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which  is  the  convolution  of  r\e  ~ivf  [rye  “*)  with  JoCPoTo* v)  and  can  be  implemented  with  an 
FFT. 

The  strength  of  this  technique  is  the  efficiency  with  which  the  Fourier  transform  can  be 
implemented.  Its  weakness  stems  from  the  requirement  that  /(r)  be  sampled  evenly  in  e~*.  In 
order  to  obtain  the  «*mplwig  density  necessary  to  represent  a  function  near  the  origin  it  is  likely 
that  such  a  density  of  points  is  necessary  to  represent  the  function  at  larger  ranges  that  the  com¬ 
putational  savings  are  lent.  Also,  the  presence  of  the  gain  factor  e  ~2v  might  be  a  severe  problem 
for  the  region  0<  r  <  1.  It  is  unlikely  diet  this  transform  technique  would  work  efficiently  for 
functions  with  more  than  moderate  range-bandwidth  products. 

In  the  next  section  we  present  another  Hankel  transform  algorithm  that  exploits  the  compu¬ 
tational  efficiency  of  the  FFt  through  a  change  of  coordinate  system.  It  requires  only  that  the 
function  be  represented  on  a  square  root  grid,  however. 

III.  7)  The  Projection-Slice  Method 
a)  Overview 

In  Section  II. 2  we  related  the  Hankel  transform  to  the  two  dimensional  Fourier  transform 
of  a  circularly  symmetric  function.  We  showed  that  the  Hankel  transform  can  be  obtained  by 
first  farming  the  projection,  or  Abel  transform: 

P(r)  -  A  f(r)  -  2f/(V^p)dy  -  2f  (1) 

o  kl  V£2-r2 

which  is  then  followed  by  a  Fourier  transform: 

9  9 

F(p)  -  Jf(r)Jfor)rdr  -  S-Jp(r)e^dr  (2) 

o  -» 

Oppenhdm,  Frisk,  and  Martinez  [11  ]  suggested  that  the  computational  efficiency  of  the  FFT  be 
exploited  by  implementing  the  Hankel  transform  as  a  numerical  projection  followed  by  an  FFT. 
Previously,  however,  the  projections  were  expensive  to  compute,  requiring  on  the  order  of  n2 
operations  and  function  evaluations  or  interpolations.  As  part  of  this  thesis  a  computationally 
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effirient  method  for  computing  the  projections  (the  Abel  transform)  was  developed.  [12  ]  When 
followed  by  an  FFT  the  result  is  a  computationally  efficient  Hankel  transform  requiring  on  thf 
order  of  N*logN  operations. 


b)  The  Abel  transform 

We  consider  the  Abel  transform  shown  below: 


P(r)~A-f(r)m  2 


i\  Vf-r1 


As  suggested  by  Bracewell  [1,2],  we  write  this  transform  as  a  convolution  by  defining: 


h(r)» 


0  r&0 

1  r<0 

—r 


and 


p(r)  »p(v7)  r*0 

/»  «/(V7)  ra:0 

which  leads  to  the  convolution  formula: 

*(')-/(')•  h(r) 

The  Abel  transform  of  /  (r  ),  p  (r),  is  determined  by: 

P(r)  =P(r*) 

Bracewell  [2]  proposed  evaluating  (6)  by  shifts  and  sums. 

Because  the  Fourier  transform  of  k  (r  )  is  the  known  analytic  function: 

H  (v)  =  Iji  for  all  v 

p(r)  can  be  determined  in  principle  by  means  of  the  Fourier  transform: 


(1) 


(2) 


(3) 

(4) 

(5) 

(6) 


(7) 


P(r)  -  fP(v)l±L±,i2».dv 

—m  *  vV 


(«) 


where  F(v)  is  die  FT  of  /  (r).  Unfortunately  the  singularity  at  v  «  0  makes  this  function  dif- 
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ficult  to  represent  in  the  computer  and  is  responsible  for  the  long  tails  of  p(r)  which  cause  alias¬ 
ing  problems  when  the  convolution  is  implemented  with  a  Fourier  transform.  We  avoid  these 

difficulties  by  removing  the  factor  from  the  numerical  part  of  the  transform  (8)  in  such  a 

way  that  the  remaining  function  is  as  well  behaved  as  F  (v)  within  the  numerical  portion  of  the 
transform.  The  singular  part  of  the  transform  is  performed  analytically  and  added  in  after  die 
numerical  processing. 

To  this  end  the  transform  (8)  is  written 


P(r)  -  r^el2^dv  + 


1 ±i  1  ./2WII 


(9a) 


!±i 

2 


b.v.il~yh +  /(Q)J ,-(>-««>  rfv  +(9b) 

+  /  I2««  J  v  +  F  (0)  f  (» V 


Where  b  is  a  parameter  analogous  to  the  real  part  of  the  exponential  in  a  Laplace  transform. 
Our  choice  of  b  is  described  later. 


The  integrands  in  the  first  and  third  integrals  of  F.quation  (9h)  do  not  have  singularities  at 
Vs1  0.  Because  both  the  numerator  and  denominator  of  these  integrands  approach  0  as  v 
approaches  0,  they  can  be  evaluated  by  1’  Hospital’s  rule  to  show  that  as  v  approaches  0  they 
approach  F(0). 

Upon  defining 


f/(0)  v-0 

{V)  "  l  (F (v)  -F(0)  e -* l*’i(l-V/v))/V/v  otherwise 
and  performing  analytically  the  two  integrals  that  do  not  depend  on  F  (v)  we  have 


(10) 


P(')' 


1+/ 

T 


dv  +  /(0) 


_ _ 1 _ _  1 

2w/r  b-2itir  V*  +2ir ir  4  +  2ir<r 


(ID 


L  (v)  was  defined  in  (10)  such  that  L  (0) (0)  and  L  (v)~  F  (v)  for  large  v.  The  parame¬ 
ter,  b  ,  was  chosen  so  that  L  (v)  moves  smoothly  between  its  limits.  If  b  were  set  to  0,  L  (v) 
would  have  a  DC  term  that  would  transform  to  an  impulse.  Theoretically  this  would  be  canceled 
by  the  singularities  in  the  portion  of  the  transform  performed  analytically  (see  equation  above) 
but  computational  errors  would  certainly  cause  problems.  If  b  were  infinite,  L  (v)  would  suffer 
from  the  lA^v  singularity  at  die  origin,  b  is  chosen  to  smooth  out  die  singularity  somewhat 
between  these  limits.  We  have  been  using  values  of  b  such  that  e~hr  has  decayed  to  e  after 
roughly  6  samples  of  F  (v). 

We  present  three  examples  of  functions  processed  with  the  Abel  transform  algorithm 
described  above. 


Example  1  (a) 


<  1 
=  1 
>  1 


(12) 


The  first  example  is  the  transform  of  the  pillbox  function 

fir)  *  {.5  r 
vO  I  r 

1024  samples  of  this  function  were  generated  on  a  grid  with  T=  1/32  and 

transformed.  The  result  is  shown  in  Figure  m.7b.l  as  dots  superimposed  over  a  solid  line  which 
is  the  transform  computed  analytically  (  2^1-r2  ).  The  output  matches  the  analytic  solution 
well. 


Example  2  (a) 

For  the  second  example  we  transformed  the  function 

^p-»ir)  (13) 

where  w  {r  )  is  a  Hanning  window.  2048  samples  on  a  VnT  grid  with  T-  1/2  were  input  (Figure 


HI. 7b. 2).  Figure  m.7b.3  diows  die  output  (dots)  superimposed  over  the  correct  transform  (solid 
line).  The  correct  Abel  transform  was  computed  by  evaluating  the  Fourier-Bessel  series  [11]  to 


Figure  III. 7b.  1  Superposition  of  the  numerically  generated  Abel  transform  of  a  pillbox  (dots) 
and  the  correct  transform 


Figure  IJL7h.2  Hi*  input  W {nT)  with  F-y  and  *  *  0, 1,  2,  •  •  •  ,  2047  where 

W(nT)  is  a  Hanning  window. 


obtain  a  slow  but  accurate  Hank  el  transform  of  the  windowed  input.  The  Hankel  transform  was 
then  inverse  Fourier  transformed  to  generate  the  Abel  transform.  In  the  absence  of  the  window 
the  result  would  have  been  sin(r)/r. 


The  output  is  coincident  with  the  correct  solution. 

Example  3  (a) 

For  the  third  example  we  again  transformed  2048  points  of 

^p-w(r)  (20) 

on  the  grid  ~^nT  but  now  T  was  chosen  to  be  4.  This  input  is  shown  in  Figure  m.7b.4. 
Increasing  the  sampling  interval  reduced  the  effect  of  windowing  on  the  input  because  a  greater 
range  of  the  function  was  represented  but  it  also  increased  the  sampling  interval  on  the  output. 
Figure  m.7b.5  shows  the  output  (dots)  superimposed  over  the  correct  transform.  Again,  there  is 
no  discernible  error. 


c)  The  Hankel  Transform 

To  complete  the  Hankel  transform  it  is  necessary  to  Fourier  transform  the  projection 
obtained  from  the  Abel  transformer.  Unfortunately  this  is  available  on  a  y^nT  grid  and  not  the 
even  grid  required  by  the  FFT.  To  generate  p  (r)  on  an  even  grid  it  is  necessary  to  interpolate. 
If  a  simple  interpolation  scheme  is  used,  like  sample  and  hold  or  Linear  interpolation,  the  result 
will  be  generated  rapidly  but  may  suffer  some  degradation.  If  a  more  sophisticated  interpolator 
is  used,  better  results  can  be  expected  but  at  the  expense  of  greater  computation  time.  Because 
the  interpolation  is  from  an  uneven  grid  to  an  even  grid  (and  not  die  reverse)  it  is  difficult  to 
characterize  die  error  theoretically  beyond  the  fact  that  the  finer  the  initial  grid  the  better  die 
results.  We  complete  the  Hankel  transform  of  the  examples  presented  above  using  linear  interpo¬ 
lation  to  generate  the  uniform  grid. 
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Example  1  (b) 

Figure  m.7c.l  shows  the  result  of  using  a  FFT  on  the  linearly  interpolated  projection  generated 
in  Example  1  (a).  The  dots  are  the  calculated  output  and  are  superimposed  over  the  analytic 
solution  (solid  line).  The  agreement  is  excellent.  The  time  required  to  perform  the  total  Hank  el 
transform  (1024  input  points  to  1024  output  points),  including  the  required  linear  interpolation, 
was  less  than  31  seconds  on  a  PDF  11-55  with  a  floating  point  processor. 

Example  2  (b) 

Figure  m.7c.2  shows  the  result  of  Fourier  transforming  the  linearly  interpolated  output  of 
Example  2  (a).  Again  the  dots  represent  the  output  of  the  Abel-Fourier  scheme  and  the  solid 
line  is  the  Hank  el  transform  as  computed  by  the  Fourier-Bessel  series.  [11]  The  agreement  is 
excellent. 

Example  3  (b) 

Figure  HI.  7c.  3  compares  the  result  of  Fourier  Transforming  the  linearly  interpolated  output 
of  Example  3  (a)  (dots)  with  the  correct  transform  (solid).  Significant  distortion  is  apparent  in 
this  transform.  Since  the  output  of  the  Abel  transform  in  example  3  (a)  essentially  equals  the 
output  in  example  2  (a)  (the  correct  projection)  except  for  the  sampling  interval,  we  can  associ¬ 
ate  this  distortion  with  the  linear  interpolation  performed  before  the  FFT. 

d)  Discussion 

We  have  found,  as  indicated  by  the  examples  above,  that  the  Abel  transformer  works  well. 
When  its  output  can  be  successfully  interpolated  and  is  followed  by  a  FFT  the  result  is  a  fast, 
accurate  Hankel  Transform  as  illustrated  by  examples  1  and  2.  As  the  spacing  between  output 
samples  of  the  Abel  transformer  is  increased,  the  suitability  of  a  simple  interpolation  scheme 
becomes  suspect.  Example  3  was  chosen  to  illustrate  the  effect  of  inappropriate  interpolation  on 
the  resulting  Hankel  transform.  At  this  point  it  would  be  prudent  to  determine  the  validity  of  a 


Hank  cl  transform  performed  with  this  algorithm  by  comparing  the  output  for  inputs  of  different 
grid  spacings. 


e)  Summary 

The  procedure  for  performing  the  Hankel  transform  H-f(r)  =  FH(v)  is  summarized 
below. 

1)  generate  /(r  )  *  / (V7 ) 

2)  Fourier  transform  to  obtain  F  (v) 

3)  generate  I(v)  *  F(v)  -  #  (0)«  1  * 1  (1  -  V^)/V^ 


4)  perform  the  inverse  Fourier  transform  and  add  in  the  analytic  terms: 


P(r) 


1FT  -L  (v)  +  F(Q) 


r 


_ _ 1  _  i^rr _ 1  1) 

ilvr  b  -2ir ir  ^b+l-ttir  b  +2itir  J J 


5)  interpolate  to  an  even  grid  p  (r  )  =  p  (r2) 


6)  Fourier  transform  to  obtain  the  Hankel  transform 

Fg(v)  =  FT  -p  (r  ) 

Each  of  steps  2)  thru  6)  requires  no  more  than  the  order  of  N  log(N)  operations.  There¬ 
fore  the  total  transform  can  be  accomplished  on  the  order  of  N  log(N)  operations.  Direct  compu¬ 
tation  of  projections  from  the  2  dimensional  circularly  symmetric  function  would  have  required 
at  least  N  function  evaluations  and  N  sums  for  each  of  N  points  before  the  final  FFT,  which 

leads  to  an  algorithm  requiring  on  the  order  of  N2  operations.  Therefore  for  sufficiently  large  N 

•  « 

tins  method  of  eaimiariwg  the  projections  can  provide  a  considerable  advantage  in  qpeed. 

Steps  1)  and  S)  above  indicate  that  in  two  places  interpolations  may  be  required.  In  many 
cases,  however,  the  function  to  be  transformed  can  be  generated  initially  on  a  grid  evenly  spaced 
inV7.  Further,  /(VnT)  ,  as  required  by  this  algorithm,  has  the  desirable  feature  that  equal 
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areas  of  the  underlying  two  dimensional  function  /  (r  ,8)  are  represented  between  samples,  if 
stationary  white  gaussian  noise  (SWGN)  corrupts  the  measurement  of  /  (  V nT  )  n  =  0,1,  •  ■  • 
then  the  Hank  el  transform  on  a  grid  will  be  corrupted  by  SWGN  (corruption  of  equal  areas 
of  the  underlying  two  dimensional  function  produces  corruption  of  equal  areas  of  the  underlying 
2-D  FT).  This  is  not  true  if /(nT)  is  corrupted  by  SWGN. 

To  implement  a  Hank  el  transform  using  this  method  it  is  necessary  to  perform  the  interpo¬ 
lations  of  step  S).  Because  of  die  speed  of  the  Abel  transform  portion  of  algorithm  we  have 
found  it  sufficient  in  many  cases  to  simply  generate  an  over  sampled  version  ofp(r)  and  to  use 
linear  interpolation  to  obtain  p  (r  )  . 

f)  Conclusion 

For  many  applications  this  method  of  calculating  the  Abel  transform  appears  to  permit  the 
efficient  calculation  of  the  Hankel  transform  for  large  data  files.  This  transform  method  is  par¬ 
ticularly  appropriate  for  the  evaluation  of  the  Sommerfeld  integral,  in  which  the  oscillations  of 
the  kernel  increase  with  the  independent  variable.  As  a  general  transform  method  issues  of 
representation  on  a  ^r"  grid  must  be  further  explored.  Because  of  the  equal  area  property 
described  earlier  for  f  )  and  because  die  speed  of  this  algorithm  permits  oversampling  in 
P  (r*)  it  is  not  expected  that  these  issues  will  pose  serious  problems.  It  appears  that  the  V/T  g^d 
is  of  fundamental  importance  in  the  Hankel  transform. 

DI.8)  Summary 

In  this  chapter  we  presented  a  number  of  numerical  techniques  for  evaluating  the  Hankel 
transform.  No  one  technique  is  ideal  for  all  situations.  When  the  value  of  the  transform  is 
desired  at  only  a  few  points,  the  Fourier-Bessel  series  or  the  backsmear  method  is  appropriate.  If 
?**d  is  extremely  important,  acoiracy  is  not,  and  the  transform  is  not  needed  for  «m«n  argu¬ 
ments,  then  the  asymptotic  method  is  justified.  If  the  input  function  and  its  Abel  transform  can 
be  well  represented  on  a  square  root  grid  (which  is  the  case  for  functions  which  increase  in 


C 
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complexity  with  range)  then  the  Hank  el- Abel  (or  projection-dice)  method  is  a  good  choice. 

Having  established  the  properties  of  the  Hankel  transform  and  examined  its  numerical 

implementation  we  are  now  ready  to  consider  uang  it  to  generate  synthetic  data. 
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CHAPTER  IV: 

SYNTHETIC  DATA  GENERATION 


IV.l)  Overview 


The  generation  of  high  quality  synthetic  pressure  fields  is  an  important  branch  of  acoustic 
research.  Because  present  methods  can  only  approximately  compute  the  fields  associated  with  a 
point  source  for  complex  environments,  simplified  environments  are  often  considered  for  which 
fewer  approximations  must  be  made.  One  important  environment  which  lends  itself  well  to 
analysis  but  which  has  sufficient  complexity  to  be  of  interest  for  practical  problems,  is  the  horizon¬ 
tally  stratified  environment.  It  is  an  excellent  model  for  the  conditions  present  in  the  deep  ocean 
over  an  abyssal  plain,  and  consequently  of  direct  interest  to  us.  Currently,  techniques  for  calcu¬ 
lating  synthetic  fields  arising  from  a  CW  point  source  in  this  environment  exist  in  the  literature 
[7,11,3  1  These  techniques  are  based  upon  the  numerical  evaluation  of  the  Sommerfeld  integral  [1 
],  for  which  two  major  computational  efforts  are  required.  First,  the  plane  wave  reflection 
coefficient  for  the  bottom  profile  must  be  numerically  generated.  For  this  the  propagator  matrix 
method  [12,8  ]  is  used.  Second,  the  pressure  field  is  computed  as  the  Hankel  transform  of  the 
depth-dependent  Green’s  function  (which  is  simply  derived  from  the  plane  wave  reflection 
coefficient).  Typically,  in  these  techniques  many  of  the  degradations  associated  with  the  numeri¬ 
cal  evaluation  of  the  Hankel  transform  are  not  carefully  controlled.  In  this  chapter  we  exploit  the 
properties  of  the  Hankel  transform  derived  in  Chapter  II  to  carefully  control  these  errors.  We 
will  show  in  Section  (IV  3a)  that  a  major  source  of  numerical  error  is  aliasing,  which  becomes 


important  because  asymptotically  the  fields  decay  only  as  —  }  We  associate  this  slow  rate  of  decay 

1 


with  the  source  singularity, 


*,  in  the  depth-dependent  Green's  function  and  show  how 


r*  VSF*?’ 

to  separate  this  portion  of  the  computation  from  the  numerically  computed  Hankel  transform. 

t)  Asmmias  tot  the  present  that  t here  are  no  trapped  modes  associated  with  low  tpaod  layers  within  the 
bottom. 


1 


The  remaining  numerical  calculations  are  significantly  iess  degraded  by  aliasing  and  are  well 
behaved  in  general.  They  remain  in  the  form  of  a  Hankel  transform  for  which  we  can  exploit  the 
computation  efficiencies  now  available  [9  ]  The  result  is  a  new  hybrid  procedure  which  is  compu- 
tationally  efficient  and  significantly  more  accurate  than  existing  methods.  This  hybrid  scheme  is 
illustrated  with  examples  of  synthetically  generated  fields. 


In  Section  (IV .3b)  we  discuss  the  difficulties  associated  with  numerically  evaluating  the 
Hankel  transform  of  the  depth-dependent  Green's  function  when  slow  speed  layers  are  present  in 

the  bottom  which  give  rise  to  proper  modes.  Proper  modes  are  associated  with  the  —r^ . - 

“ifj 

singularities  (with  near  the  real  axis)  in  the  depth-dependent  Green's  function  and  contribute 


terms  to  the  field  which  decay  asymptotically  as  — yr.  This  very  slow  decay  in  the  field  causes 


severe  aliasing  problems  when  it  is  calculated  using  a  numerical  Hankel  transform  algorithm.  In 
Section  (IV Jb)  we  show  how  to  separate  the  effects  of  proper  modes  from  the  numerical  calcula¬ 
tions  by  performing  part  of  the  Hankel  transform  analytically.  We  make  this  separation  in  such  a 
manner  that  the  portion  of  the  field  assigned  to  the  analytical  calculations  is  exact  and  finite  for 
all  ranges.  This  makes  it  possible  to  numerically  calculate  the  residual  numerical  contribution  to 
the  field  accurately  and  add  the  result  to  the  analytically  determined  expression.  The  result  of  the 
total  hybrid  algorithm  is  a  field  which  is  accurate  for  all  ranges  and  which  can  accurately  include 
the  rffects  due  to  proper  modes  arising  in  the  presence  of  slow  speed  layers.  We  present  an  exam¬ 
ple  of  a  field  generated  synthetically  with  this  total  hybrid  method  and  show  how  the  result  is 
significantly  better  than  what  would  have  been  achieved  without  removing  the  effect  of  the  poles. 

In  this  chapter  we  also  develop  a  numerical  implementation  of  the  propagator  matrix  method 


for  generating  the  plane  wave  reflection  coefficient  that  is  well  behaved  numerically.  We  begin 
the  chapter  by  describing  the  computation  of  the  plane  wave  reflection  coefficient  by  means  of  the 
Thomson-HaskeU  method  [12,8  ]. 


IV .2)  The  Propagator  Matrix  Approach  to  Generating  the  Plane  Wave  Reflection  Coefficient 


a)  The  Method  in  Principle 
I)  Overview 

To  calculate  the  plane  wave  reflection  coefficient  we  coniider  the  responte  of  a  layered  bot¬ 
tom  to  an  incident  plane  wave  aa  shown  in  Figure  IV2ai.l.  Within  the  n1*  isovelocity  layer  we 
express  the  field  as  the  vector: 

P(z) 

U(z)  a) 

where  Pn(z)  is  the  pressure  in  the  n,h  layer  and  Un  (r)  is  the  normal  component  of  the  velocity. 
We  have  chosen  this  representation  because  P  (z )  and  U(z)  are  continuous  in  r ,  even  across 
layer  interfaces.  In  the  discussion  which  follows  we  will  suppress  the  time  and  radial  dependence 
of  the  field  because  they  are  the  same  in  all  layers. 

In  the  propagator  matrix  approach,  the  impedance  at  the  bottom  layer: 

P  (zy  +l) 

is  available  from  the  material  parameters.  In  principle  this  impedance  is  used  to  determine  the 
reflection  coefficient  at  the  top  interface  in  three  steps.  First  the  field  at  the  top  interface  is 
related  to  the  field  at  the  bottom  interface  by  the  propagator  matrix: 

>(*<>)]  [F(zy)l 

f/(z0)J  "  *  [u(zN)\  (3) 

Next  the  incident  and  reflected  pressure  waves  at  the  surface  are  related  to  the  field  at  the  top 
interface  and  then  to  those  at  the  bottom  through: 

f(*o)l  pH**)]  k  an]  pH**)] 

P-fl\  "All/(*o)j  “A *>(**)]  "  [«ai  «»]  [U(zN)\ 

Finally,  the  reflection  coefficient  is  calculated  in  terms  of  the  impedance,  {y +l,  using 


(4) 
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so  that 


P  -.0  _  *11  +  Ctf+la12 


(6) 


Xo  “21  +  Cy+la22 

U)  Th*  Propagator  Matrix 

The  essential  element  of  this  approach  is  the  propagator  matrix,  $,  of  Equation  IV2a.l3.  In 
this  section  we  review  its  derivation. 

Within  any  isovelocity  layer,  the  held  can  be  considered  as  the  superposition  of  a  positive 
and  a  negative  traveling  wave.  The  pressure  field  is  given  by 


P(z )  «  ?+«*•*  + 


(1) 


The  normal  component  of  velocity,  V  (2 ),  is  related  to  P  (z )  through  the  telegraph  equations  [2  ] 

For  the  non-normal  case  we  use  — —  «  “  -iwpf/  which  implies  that: 

dz  it 


U(z)~-r —P+t 

tup  tup 


or  defining  Y  q  *  - : 

up 


U (s)  ■  Y^P**'*'*  -YoP-e'*’’ 
In  matrix  form  Equations  (1)  and  (3)  become: 


(2) 


(3) 


P(z) 

X 

U(z) 

lync*4’* 

P_ 

(4) 

X*  1) 

\ 

F(rj) 

If 

utod. 

is  known  at  tome  point  in  the  layer  then 

?<*2> 

can  be  computed  in  principle  by 

X 

1 

1 

>(*1) 

>(*2) 

X 

inverting  Equation  (4)  to  find 

P. 

»  < 

in  terms  of 

and  then  calculating 

X*2). 

from 

/>_  • 

Combining  these  operations  into  one  step  gives: 


>(*a> 

2) 


A»s 


y0si*,,‘  -y 


-1 

>(*1) 

u(tx) 

(5) 


which  gives  us: 
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X*z) 

V(*a) 


cosAj  (Z2~Z  l)  -^-sin*,(z:-zi) 

'  0 


X*t) 

U{*  i) 


(6) 


posin*«(*a“*i)  «**»(*  r"*i) 

The  values  of  kt  and  Y  0  are  functions  of  the  material  parameters  of  the  layer  under  considera¬ 
tion.  In  particular  if  c*  is  the  sound  speed  in  layer  n,  p*  its  density,  kr  the  horizontal  wave 
number  of  the  incident  plane  wave  (by  Snell’s  law  common  to  all  layers),  and  u  the  temporal  fre- 

Jfc 

quency  of  the  CW  source,  then  k„  ■  — ,  k,  »  \/k~—k~t  and  Yq  m  — — . 

c„  <up 

To  indicate  explicitly  the  dependence  on  the  material  properties  of  the  n,k  layer  we  write: 


X*2> 

U(Z2) 

when  z  j  and  z  i  are  both  within  layer  n . 


X*  i) 


(7) 


To  calculate  the  field  at  the  top  interface  in  terms  of  the  field  at  the  bottom  interface,  as 
shown  in  Figure  IVlai.l,  we  can  use  the  previous  discussion  which  was  applicable  only  to  a  single 
layer,  to  relate  the  field  at  z,  to  the  field  at  z,_1: 


X*,-i) 

fffes-t). 


X*.)‘ 


>(*,) 

UizK) 


We  then  iterate  the  procedure  through  all  the  layers  to  find 


(8) 


X*o) 
P  (*o) 

b)  Numerical  Implementation 


X*„) 


■ 


(9) 


I)  TAe  modified  propagation  matrix 

The  bulk  of  computation  associated  with  the  propagator  matrix  approach  is  the  accumulation 
of  the  matrices.  ^1^2  1  • '  When  these  are  accumulated  on  the  computer,  the  actual  opera¬ 
tion  is  •••$// ))).  It  is  possible  for  the  scale  of  the  accumulated  product  to  differ 

dramatically  from  any  particular  .  Because  of  the  limited  dynamic  range  in  the  computer  it  is 
advisable  to  scale  terms  to  make  them  comparable  before  accumulation.  Fortunately  the  final  cal- 
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culation  for  the  reflection  coefficient  depend*  only  on  ratio*  of  element*  in  4>.  For  this  reason  we 
normalize  each  of  the  to  that  its  largest  value  equals  1.  This  procedure  alone  could  cause 
another  problem  stemming  from  the  different  scales  in  general  for  P  and  U ,  which  is  due  to  their 
different  units.  To  bring  P  and  U  into  the  same  units  we  do  not  actually  relate 

[/>(*„) 

t0  [t/M  (1) 

but  rather  consolidate  units  by  multiplying  the  normal  component  of  velocity  by  the  characteristic 


impedance  of  that  layer.  Therefore  we  actually  calculate: 


F(*0) 
iou (*o) 


aM  bH  1  \PM 
lNbN  U  (2W  ) 


where 


.  _  «»P i 
T 


a,  *  cos 

bi  ■  -isin  kijizi-ti- 1) 

U)  Relation  of  the  modified  propagation  matrix  to  th*  incident  and  reflected  waves 

We  now  relate  the  field  variables  to  the  incident  plane  wave  and  the  resulting  reflected  plane 
wave  by  slightly  modifying  Equation  (lVJ2ali.4).  We  assume  that  the  top  interface  is  at  z  *  0  so 


>(0)  fl  0]  [l  1  1  \P+i 

:Co^(0)J  “  10  Co]  [Yo  -Y  oj  [/»-; 


1  ± 


l  yo  \l  0] 

“  2  i  dL  [0  Y0\  [| 


P(0) 
toU  (0) 
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By  defining: 

a,  bt 
ki*>i  iiOi 

and  using  Equation  (IV  2b  .i  2)  we  have: 


$U  $12 

V 

,$21  $22. 

-  n 

im  1 

(3) 


(4) 


p+>o  i  fl  ll  f$u  $12]  p*(2iv) 

mp -a  2  L1  ”lJ  [$21  $22]  [Cw u  (** ), 

We  now  need  to  use  the  fact  that  the  pressure  and  velocity  fields  in  the  last  layer  are  made  up  of 
only  positive  traveling  waves  so  that  (referring  to  Equation  IV2aii.4)  PN  +l  -  P  (zN )  and 

tfjv+i  ■  “—/*  (*y  )  we  have 

t-l 


If  we  now  use 


/*+.<> 

_  1 

$U+$21  $tt+$j2 

I 

'-.e. 

2 

,$11“$21  $12“$22. 

Cat+i 

(5) 


(6) 

VW  -t-l 

we  have  the  reflection  coefficient 

T  «  P  ",0  -  +  t/f  +i($i2~$22) 

P+,0  $u‘*'$2i  +  ijv+i($:2+$22)  ' 

Equations  (IV2bii2)  and  (7)  show  that  this  approach  uses  only  the  ratios  of  the  impedances 
in  adjacent  layers  and  never  the  impedances  themselves,  "uese  ratios  are  much  better  behaved  in 
general  than  the  individual  impedances.  For  this  reason,  because  the  use  of  P  and  Y^U  instead  of 
P  and  U ,  and  because  of  the  scaling  of  the  layer  propagation  matrices  this  implementation  of  the 
propagator  matrix  approach  has  good  numerical  properties. 


c)  Selected  Properties  of  the  Rtf  faction  Coefficient 

la  Figure  IV2c2  we  present  a  perspective  plot  of  the  log  magnitude  of  the  reflection 
coefficient  as  a  function  of  kr  for  the  bottom  of  Figure  IV2c.l  calculated  using  the  numerical 
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f  *  220  Hz 
|z  +  z0 1  a  2 m 
k0s  .8975979  ftT1 


C0  8  1540  m/s 
p0 *  1  g/cm3 

/  /  /  /  .  /  ,•  . 

C,  -  1493.8  m/s 

yO,  =  1.5  g/cm3 

C2  *  1700  m/s 
pz  *  2.0  g/cm3 


Figure  TV.2c.l  The  bottom  parameters  used  to  generate  the  reflection  coefficient  shown  in  Fig¬ 
ure  IV.2c.2 
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algorithm  just  described.  The  reflection  coefficient  is  displayed  both  for  those  horizontal  wave 
numbers  corresponding  to  real  angles  of  incidence  (0  <  Re(Jfcr )  <  *0  where  k0  is  the  water  wave 
number)  and  for  horizontal  wave  numbers  corresponding  to  complex  angles  of  incidence 
(*o  <  ))•  The  complex  angles  correspond  to  evanescent  waves.  Single  evanescent  waves  do 

not  carry  any  time  average  power  flow  (their  Poynting  vector  is  imaginary)  and  consequently  the 
magnitude  of  the  reflection  coefficient  is  not  limited  to  be  less  than  one  in  the  evanescent  region, 

*0<*r  [«) 

In  Figure  W2c2  a  pole  is  apparent  in  the  reflection  coefficient  on  the  real  kT  axis  in  the 
evanescent  region.  This  on  axis  pole  corresponds  to  a  proper  mode  propagating  in  the  low  speed 
layer  within  the  bottom.  Other  off  axis  poles  corresponding  to  leaky  modes  are  also  apparent  in 
the  reflection  coefficient.  A  discontinuity,  or  cut,  can  be  seen  extending  from  k2  along  the  real  kr 
axis  to  infinity.  This  is  the  branch  cut  extending  from  the  branch  point  at  k2.  Another  cut  extend* 
ing  from  k0  to  infinity  falls  on  this  same  line  and  is  therefore  not  apparent.  The  origin  of  these 
branch  points  and  cuts  can  be  found  in  our  derivation  of  the  reflection  coefficient  where  we  asso¬ 
ciated 
« 

«‘v kj-v*  with  (1) 


e~i\/kj-k*i  with  p_  @) 

Clearly  the  roles  of  P  +  and  P  _  would  be  reversed  by  changing  the  choice  of  sign  for  the  square 
root.  For  incident  and  reflected  wave  this  would  correspond  to  inverting  the  reflection  coefficient 
(if  no  other  waves  were  affected).  The  two  sheets  corresponding  the  the  branch  point  at  k0 
reflea  the  two  choices  of  sign  for  the  incident  wave.  We  have  displayed  the  choice  associated  with 
positive  real  power  flow  for  the  incident  wave. 

In  the  intermediate  layers  such  as  layer  1  of  this  example,  changing  the  role  of  and  /*_ 
would  not  affect  the  reflection  coefficient.  For  the  intermediate  layers,  the  physical  problem  does 
not  name  (or  distinguish  between)  forward  and  backward  traveling  waves.  Consequently  there  are 


C 


A 


no  branch  points  associated  with  intermediate  layers. 

If  the  opposite  sign  were  chosen  for  the  square  root,  V*w+i associated  with  the 
emerging  wave,  P  from  the  bottom  of  the  stack  of  layers  (into  the  isovelodty  half  space)  the 
direction  of  energy  flow  associated  with  that  wave  would  change.  Unlike  the  intermediate  layer, 
there  is  no  returning  wave  in  the  isovelodty  half  space.  Consequently,  the  physical  problem  would 
change.  For  this  reason  we  tee  a  branch  point  at  k# +1  reflecting  the  two  different  'physical'  prote¬ 


in  Figure  lV2c2  we  have  chosen  to  display  the  Riemann  sheet  for  which  both 
Re(v£F*7)  >  0  and  Re(V^jv7i"“*7)  >  0-  On  this  sheet  only  waves  with  real  power  flow  in 
the  positive  direction  are  associated  with  P  +.  This  constrains  our  incident  waves  to  be  those  with 
power  flow  from  the  source  to  the  layered  bottom  and  specifies  that  there  is  no  power  flow  return* 
ing  from  infinity. 

When  we  perform  the  integrations  discussed  later  we  must  choose  which  side  of  the  cuts  to 
integrate. upon.  For  reasons  of  convergence  we  choose  the  side  for  which  ImCVA/-*,2)  >  0 
when  j m  0  and  N  +1,  is  satisfied.  Consequently,  whenever  we  integrate  the  reflection 
coefficient  in  the  complex  k,  plane,  we  always  satisfy  both  Re(\A/~*rz)  >  0  and 
>  0  for  j  -  0  and  N  +1. 


IV  4)  Evaluating  the  Sommarftld  Integral 

Once  the  plane  wave  reflection  coeffideat,  T  (k, ),  has  been  computed  it  is  necessary  to 
evaluate  the  Sommerfeld  integral: 


Pg(r)  -  f  -  ^^r(krWyy*F*',*'''j<trkr)krdkr  (1) 

o  V*o-V 

in  order  to  compute  the  reflected  pressure  field.  The  Sommerfeld  integral  is  in  the  form  of  a 


Hank  el  transform  of  the  depth  dependent  Green’s 


function, 


G(k,jj0)  ■  ■  r(kf )«' The  properties  of  the  Henkel  transform 


« 
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developed  in  Chapter  II.  In  Sections  (6)  and  (7)  of  that  chapter  we  discussed  the  effect  of  trun¬ 
cating  the  integration  at  some  finite  value,  and  the  effect  of  sampling.  The  truncation  was  accom¬ 
plished  by  multiplying  the  function  to  be  transformed  by  some  finite  length  window.  For  the  gen¬ 
eration  of  synthetic  data  we  find  that  windowing  of  the  Green's  function  is  not  an  important  con¬ 
sideration  in  general  because  when  x  +*o  >  0,  G  (k,  4  ^q)  decays  exponentially  in  k,  for  kr  >  k0. 
Except  when  *  +i0  is  very  small  we  can  integrate  Equation  (1)  until  the  Greens  function  is  negli¬ 
gible  and  truncate  at  that  point.  It  is  not  necessary  to  multiply  by  a  windowing  function. 

The  issue  of  sampling  and  the  associated  degradation  introduced  into  the  transform,  aliasing, 
can  be  very  much  a  problem  however. 


a)  The  Source  Singularity 

In  order  to  highlight  the  issues  associated  with  the  source  singularity,  —  ‘  and  the 

propagation  terms  in  the  Sommerfeld  integral,  we  first  consider  the  evaluation  of  Equation 
(IV  3.1)  for  a  hard  bottom  case  where  r(Jk,)  »  1. 

For  the  hard  bottom  case  the  pressure  field  is  given  by  the  known  integral: 


P*(r) 


f  ■* 
0 


ei\£FPi 


(2) 


Vr2+(S  +*o)a 

We  evaluated  this  integral  numerically  with  the  Fourier-Bessel  series  to  obtain  an  estimate  for  the 
field: 


0  <  r  <  A 


(3) 


where  for  this  example,  A  was  chosen  to  be  2000  and  !*  +r0l  ■  2  .  In  Figure  (TV3a.l)  we  com¬ 
pare  the  log  magnitude  of  the  result  (dots)  with  the  known  transform  (solid  curve).1  We  see  that 
the  magnitude  of  the  numerically  generated  field,  Pg  (r ),  oscillates  rapidly  in  contrast  with  the 


l)Tho  output  of  tbo  Pourlor-Bowol  win  has  boss  di^lsyod  to  twieo  its  npot  of  validity  to  bottot  Ulus- 
trsto  tbo  soutco  of  dofradttioo. 
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Figure  IV.3a.I  A  comparison  of  the  log  magnitude  at  the  reflected  field  generated  by  a  point 
source  over  a  hard  bottom  (solid  Use)  to  the  Add  numerically  computed  using  the  Fourier-Bemel 
aeries  (scatter)  which  is  shown  to  twice  its  presumed  region  of  validity 


i 

■i 

i 
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true  Hankei  transform.  As  we  will  now  show  these  oscillations  are  due  to  aliasing  in  the  numeri¬ 
cally  computed  Hankei  transform. 


in  Section  (11.7)  we  showed  that  the  effect  of  sampling  on  the  Hankei  transform  is  to 
approximately  produce  an  aliased  estimate  of  the  true  transform,  v7?*(r).  Since  for  this  exam¬ 
ple,  P/t(r)  decays  asymptotically  as  1/r,  vFf*  (r  )  decays  asymptotically  only  as  1  /vF.  What  we 
see  in  Figure  (IVJa.l)  is  given  approximately  by: 


0<  r  <  2A  l>j ,(r)l  « 


1 

v7 


Vr'+Cs+s  o)4 

When  r  is  much  greater  than  r  -fro,  Fx(r )  is  approximately: 


_  u0\/(M -')*+(»  ♦to)1 

VF-  =-r  - 


VCZA-rM^+Xo)4 


(4) 


1  1 

^iky  ^ii0l2A-rl]| 

Iv7 

vF  V5T7  JJ 

Since  we  are  in  the  region  r  <  2A  this  can  be  rewritten: 


(5) 


0<r  <2A 


.“o' 


e  _  -iky 


vT  ~  v£T7e 


We  can  write  Equation  (6)  in  terms  of  the  desired  transform  and  a  modulation  term  as: 


1  1 

f  l 

ikf^A 

*y  .  «{  «  r!n  t  „ 

Iv7 

IvT  vSwj 

v\ 

Which  upon  defining  «(r  ) 


i^)i 


= 


1  tfr) 

r  VT 


+  2i^p-sin  *or 


(6) 


(7) 


(8) 


When  r  «  2 A  *(r)  is  small,  so  that  the  magnitude  of  P&(r)  appears  as  roughly  the  correct 
transform  with  a  modulation  term. 


We  note  at  this  point  that  if  we  had  sampled  the  output  of  our  transform  at  an  integral  mul¬ 
tiple  of  2ir/*o  we  would  not  have  teen  these  oscillations.  At  this  sampling  rate  the  cosine  would 
have  appeared  as  a  DC  offset  in  the  magnitude  of  the  pressure  field.  If  the  output  nwpiing  rate 
were  near  but  not  exactly  an  integral  multiple  of  2w/Jk0  the  cos(k</)  would  have  appeared  as  a 
low  frequency  modulation  because  the  sampling  is  in  effect  demodulating  the  cosine  down  to  a 
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low  (but  not  bow  sere)  frequency.  Tbit  remit  it  an  important  one  because  frequeatly  pressure 
fields  are  generated  by  using  an  FFT  based  approximation  to  the  Hank  el  transform  (TV  .3.1)  and 
the  water  wave  number  is  the  maximum  wave  number  used  [5,4  ]  The  grid  resulting  from  such 
processing  u  exactly  an  integral  multiple  of  2n  /k  q.  Carrying  the  integration  to  higher  wave 
number  would  make  evident  the  modulation  in  the  answer  by  automatically  providing  the  output 
on  a  finer  grid. 

The  problem  of  aliasing  arose  because  the  field  being  computed  decayed  only  as  which 
forces  us  to  use  a  very  high  sampling  rate  to  properly  sample  the  Hankel  transform.  We  now  note 
that  this  -4r  decay  is  due  to  the  _1__-  singularity  in  the  Green’s  function.  It  is  well  known 

Vr  y/kF*? 

that  the  asymptotic,  or  far  field,  character  of  such  a  transform  is  determined  by  the  singularities 
of  the  kernel  over  the  path  of  integration  [10  ].  The  Green’s  function  which  is  transformed  in 
Equation  (IV  3.1)  was 


G  {kr  ji  ,r  o) 


iVtJ-*,1!*  i 


The'  asymptotic  character  of  the  transform,  P  (r  ),  is  dominated  by  the  singularity 


1 

The  integral 


(9) 


(10) 


e  a  m  r  i 
Vr1+«*  o  “\/k o  —kf 


Jo(krr)krdkr  (11) 

shows  us  that  this  singularity  is  in  fact  associated  with  the  1/r  decay  rate.  Physically  this  singular* 
ity  was  due  to  the  angular  spectrum  of  the  point  source.  The  1/r  deeay  associated  with  this  singu¬ 
larity  is  often  associated  with  the  point  source  by  noting  that  the  field  around  a  point  source  must 


decay  at  that  rate  in  a  manner  such  that  the  intensity,  which  decays  as  the  field  squared, 
integrated  over  any  three  dimensional  shell  enclosing  the  point  source,  would  not  be  a  function  of 


r. 
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The  source  of  slow  asymptotic  decay  we  have  isolated  suggests  a  procedure  for  reducing  the 
problem  of  aliasing.  We  remove  the  source  singularity  from  the  kernel,  numerically  transform 
what  is  left,  and  add  the  result  to  the  analytically  determined  transform  of  the  source  singularity. 
When  we  remove  the  singularity  we  must  do  so  in  a  manner  such  that  the  numerical  transform  we 
must  perform  is  well  behaved.  We  allow  for  a  general  F(k,)  but  at  this  time  assume  that  r(*r) 
has  no  singularities  along  the  real  kr  axis  with  asymptotic  contributions  to  the  field  capable  of 

dominating  those  of  the  singularity. 

To  this  end  we  rewrite  integral  (IV  3.1)  as: 


Akrr)k'dk' 


’/  r(k,)-r(*)j 

a  l  * 


,1 _ L 


vST-i? 


/n(k,r)k,dk,  +r(k)j 


(12) 

'/ o(k,  r  )k,dk,  (13) 


If  we  define: 


L(kr)m  [r(M  -  (14) 

so  that  L  (kr)  does  not  have  the  1  A/k o  ~k?  singularity  at  kr  *  kg.1  then  we  can  write  Equation 
(13)  as: 

*  it  VV  *■*■(« 

fair)  -  J L{k,)J0(krr)krdkr  +  r  (*)-£■  :  =•  (15) 

o  Vra+(*  +*o) 

Because  L  (kr)  does  not  have  this  singularity  along  the  path  of  integration  the  output  of  the 
numerical  transform  will  decay  at  a  rate  faster  than  1 lx.  The  asymptotic  1/r  decay  is  provided  by 


We  »aow  ia  the  appendix  that  if  the  impedance  and  It*  first  derivative  at  the  interface  is  finite  for 
kr  m  k  o  then  the 


Hm  L(kr) 

v'o 


«*p# 


where  Zi  is  the  impedance  of  the  bottom  at  k,  •  k*  <•  it  2w  source  frequency,  and  p«  is  the  density  of 
the  water.  For  an  isoveiocity  half  space  this  expression  reduces  to 


Which  is  finite. 


the  analytic  term  which  can  be  recognized  at  the  tpecular  reflection  when  r  it  very  targe  (glancing 
incidence).  These  observations  are  confirmed  in  the  examples  which  follow. 

In  the  following  examples  we  illustrate  the  generation  of  synthetic  pressure  fields  through 
the  hybrid  algorithm  implied  by  Equation  (15)  where  the  integral  is  performed  with  a  numerical 
Hankel  transform  algorithm  and  the  analytical  expression  it  the  result  of  integrating  the  singular* 

ity. 

1)  Hard  Bottom 

This  is  the  degenerate  example  because  for  r(kr)  constant,  the  entire  transform  is  per¬ 
formed  analytically.  The  result  of  the  analytic  transform  was  compared  to  the  direct  numerical 
transform  Figure  (IVJa.l). 

U)  Slow  bottom 

Figure  (IV.3aii.l)  shows  the  bottom  parameters  for  this  example.  Figure  (IVJaii2)  shows 
the  result  of  the  hybrid  calculation  (solid  line)  versus  a  direct  numerical  calculation.  The  improve¬ 
ment  is  dramatic.  Figure  (IVJaiiJ)  compares  the  hybrid  field  of  Figure  (!VJaii2),  with  its 
numerically  generated  component.  As  can  be  seen,  the  near  field  is  dominated  by  the  numerically 
generated  component.  As  range  increases  this  numeric  term  begins  to  suffer  from  aliasing  prob¬ 
lems  but  the  analytic  term  begins  to  dominate,  minimizing  the  effect  of  aliasing  on  the  computed 
field  at  large  ranges. 

ill)  Fast  Bottom 

Figure  (IV3aiii.l)  shows  the  parameters  of  the  fast  bottom  for  this  example.  Figure 
(IVJaiii.2)  shows  the  hybrid  calculation  versus  the  direct  numerical-  calculation.  Figure 
(IV.3a.iii3)  presents  the  hybrid  field  and  its  numeric  component.  The  improvements  are  similar  to 
the  fast  bottom  cate. 


b)  Poles  Due  to  Slow  Speed  Layers 
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f  =  220  Hz 
|  z  ♦  z0|  s  2  m 
k0*  .8975979  m”1 


C0S  1540  m/s 
y00=  1.0  g/cm3 

TTT/T////T//7/T 

C,  =1700  m/s 
py  -  2.0  g/cm3 


Figure  IV.3s.ULl  Parameters  of  bottom  used  for  fast  bottom  example 
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Figure  (IV3b.l)  show*  the  parameter*  of  a  slow  speed  layer  between  two  isoveiocity  half 
spaces.  These  axe  the  same  parameters  used  to  generate  the  perspective  plot  of  the  reflection 
coefficient  presented  in  Figure  (IV  2c2).  Figure*  (IV.3b.2a)  and  (IV 3b 3b)  show  the  magnitude 
and  phase  of  the  reflection  coefficient  for  this  bottom  as  a  function  of  horizontal  wave  number. 
We  see  that  for  this  example  the  reflection  coefficient  has  a  singular  point  beyond  the  water  wave 
number,  That  singularity  is  a  simple  pole  associated  with  a  proper  mode  excited  in  the  low 
speed  layer.  Such  a  proper  mode  can  appear  only  for  kg  <  k,  <  kN+  In  this  region  conserve* 
tion  of  energy  is  not  violated  because  the  waves  are  evanescent.  Proper  modes  are  generated 
when  the  low  speed  layer  acts  like  a  dielectric  waveguide.  When  this  happens  energy  diffuses  (tun* 
nets)  from  the  point  source  to  the  low  speed  layer  but  does  not  otherwise  propagate  vertically. 
Energy  from  the  field  is  now  constrained  to  decay  in  only  two  dimensions  instead  of  three  and  we 

expect  that  the  field  associated  with  the  pole  will  decay  asymptotically  as  — 4r  so  that  the  integral 

Vr 

of  the  flux  over  any  two-dimensional  ring  surrounding  the  source  remains  constant. 

Poles  such  as  this  disrupt  the  asymptotic  character  of  the  field  derived  in  the  previous  sec¬ 
tion.  As  before  we  would  like  to  analytically  determine  the  contribution  of  these  poles  and 

remove  them  as  we  removed  the  ■  1  ■  »  singularity.  To  do  so  it  is  necessary  to  evaluate  the 

integral: 


l(r  J+Z0*rj  *  / 


i  vPhE7ii+*0i 


vP=5? 


J0(k,r)krdkr 


In  Appendix  (I)  we  show  that  for  lm(kr  )  %  0  (associated  with  no  return  from  r 
/(r^+zo^,)i«  given  by: 

. -"a  -  l.-*1'-' 


(1) 
0  ) 


(2) 


where 


(3) 


|z  *  z0 1  a  2  m 
k0*. 8975979  m*1 


C0  *  1540  m/s 
p0 *  1  g/cm3 


///•it  /  /  //////.  /  /  ' 

C,  *  1493.8  m/s  a 
yo,  *  1.5  g/cm3  v 


10  m 


C2  *  1700  m/s 
pz  «  2.0  g/cm3 


Figure  IV.3b.l  Parameters  of  bottom  used  for  dow  speed  layer  example 
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The  first  integral  is  easily  evaluated  on  the  computer.  In  addition  as  r  becomes  large  the  first 
term  rapidly  approaches: 

♦»<!)* 

2?  Vr’+C* +T0?  W 

The  second  term  decays  as  1/vF  when  kfi  is  real. 

Equation  (2)  is  also  correct  tor  Im(hf  )  >  0,  but  when  Im(k,()  »  0  the  poles  no  longer 
contribute  asymptotically  as  1  /v7  because  the  Hankel  function  decays  exponentially.  Under  these 
conditions 


C(*f  )  -t,  f  ikf  -r 

- - ——e  1  t  4 

vT 


(5) 


(6) 


As  Im(Jkr  )  becomes  large  the  exponential  decay  dominates  the  1/Vr*  decay  even  over  the  finite 
range  that  concerns  us.  It  is  for  this  reason  that  we  consider  only  those  poles  near  the  real  axis 
(close  to  the  path  of  integration)  and  leave  the  others  to  the  numerical  pan  of  the  transform. 

With  /  (r  ,z  ^Zq-Jc,  )  so  defined,  the  reflected  pressure  field  can  be  written: 


MO-J 


i 


r(*'  *  2_>  2 
I  Kr  *f 


#tvP=?.»«0 'j^kr,)k'dk'  +  2c/(r  ^+*o;  kr.)  (7) 


o  y/k$-k; 

Where  the  expression  in  brackets  no  longer  has  any  poles  near  the  line  of  integration  and  so  can 
be  evaluated  as  before. 


In  order  to  remove  the  poles  as  required  in  Equation  (7)  it  is  necessary  to  determine  with 
precision  the  pole  locations,  kr>,  and  their  scales,  (a _j)( .  The  pole  locations  can  be  found  using 
standard  complex  root  finding  techniques,  though  care  mutt  be  taken  to  provide  the  root  finding 
algorithm  with  values  of  the  reflection  coefficient  on  the  Riemmaa  surface  so  that  it  appears  ana* 
lytic  except  at  isolated  singularities.  This  means  that  the  branches  chosen  for  the  square  roots 
must  be  taken  in  such  a  manner  that  a  branch  cut  is  never  placed  between  points  simultaneously 
considered  by  the  root  finder.  Once  the  root  locations  are  known,  the  scale  factors  can  be  found 
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(or  those  singularities  far  from  any  other*  by  determining  a  least  square*  fit  to: 


r(*.) 


(a-i)i 


j  -U.  N 


(8) 


V  k  2-k  2 

provided  that  the  kfj  are  taken  sufficiently  dose  to  kr.  that  T(kr)  is  well  approximated  by  just  one 
pole  in  that  region. 

If  many  poles  are  clustered  together,  they  can  be  determined  simultaneously  by  solving: 


(9) 


for  N  sufficiently  large.  If  a  pole  is  near  a  branch  cut  then  the  poles  on  the  other  side  of  the  cut, 
on  the  opposite  sheet,  and  near  the  cut  must  alto  be  considered  to  be  near  that  pole. 

Figures  (IV  ,3b  3a)  and  (IV  3b  3b)  show  the  magnitude  and  phase  of  the  reflection  coefficient 
of  Figure  (IV 3b 3a)  minus  the  pole  contribution: 

a-l 


r(*')  * 

For  this  example  o.s  -  1.689712*  10*1  i  5.027826*10" 

kr  -  9.069830* 10"1  +  i  2.488749*  10"5. 


(10) 

and 


Note  the  difference  in  scale  between  Figures  (IV 3b 3a)  and  (IV 3b 3a).  The  small  notch  visi¬ 
ble  at  Jfc,  ~  kr  is  due  to  a  small  amount  of  error  in  the  estimate  of  kr  . 

A  notable  feature  of  Figure  (TV  3b  3a)  is  the  unmasking  of  the  off  axis  seres  in  the  region 
k/t+i  <  k,  <  k0  where  previously  ir(k,)l  ■  1.  These  seros  can  be  slcarty  seen  in  the  perspec¬ 
tive  plot  of  the  total  reflection  coefficient  in  the  complex  plane  that  was  pa  eased  in  Figure 
(IV  3c  3). 

Figure  (TV3b4)  presents  the  hybrid  fleld  (solid  Uae)  verses  the  held  eel  salami  without 
removing  the  pole  from  the  reflection  coefficient  (bat  otherwise  weeing  the 
larity  as  in  the  previous  hybrid  examples).  The  spread  in  the  direedy  competed  fleU  dee  so  alien¬ 
ing  is  severe  because  aliasing  in  the  Hank  si  transform  severely  affects  a  function  that  decays  a* 
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Figure  TV.3b.3e  Magnitude  of  reflection  coefficient  for  flow  speed  layer  example  after  the  pole 
has  been  removed 
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1/VT .  The  hybrid  held  does  not  exactly  follow  the  contour  of  the  top  of  the  ipread  just  as  the 
hybrid  computations  in  the  previous  examples  did  not  exactly  follow  those  contours  when  the 
aliasing  became  severe.  Figure  (IV 3b .5)  presents  the  log-magnitude  of  the  analytically  generated 
pole  contribution  (solid  line)  and  the  remainder  of  the  held  exclusive  of  the  pole  contribution. 
The  non-pole  contribution  is  most  significant  for  short  ranges,  while  for  this  near  bottom 
geometry  the  pole  contribution  dominates  farther  out. 

The  expression  for  /  (r  ,z  +*<>&•,)  in  Equation  (2)  shows  that  the  contribution  of  the  pole  to 
the.  field  decreases  exponentially  with  lz+r0|.  In  this  example  !z+z<>l  *  2  to  emphasize  the 
near  field  behavior  associated  with  the  pole.  For  larger  values  of  !z  +r0l  the  pole  contribution 
would  be  considerable  less.  Equation  (2)  can  be  used  to  estimate  the  magnitude  of  the  pole  con¬ 
tribution  if  the  pole  location,  kr. ,  and  I  z  +z0l  are  known. 
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CHAPTER  V: 

THE  INVERSION  PROCEDURE 

V.l)  Overview 

la  Chapter  n  we  developed  die  properties  of  die  Henkel  transform.  These  properties  pro¬ 
vided  our  foundation  for  the  development  of  an  accurate  procedure  to  numerically  generate  syn¬ 
thetic  pressure  fields,  presented  in  Chapter  TV.  In  this  chapter  we  will  use  die  results  of  Chapter 
n  to  explore  the  problem  of  determining  the  plane  wave  reflection  coefficient  from  measure¬ 
ments  of  the  pressure  field  arising  in  a  horizontally  stratified  environment  in  response  to  a  CW 
point  source.  The  estimation  of  the  plane  wave  reflection  coefficient  from  measurements  of  the 
field  is  an  extremely  important  problem.  Determining  die  plane  wave  reflection  coefficient  is  an 
essential  step  in  the  inversion  of  pressure  field  data  to  obtain  the  parameters  of  the  bottom.  In 
this  context  it  is  of  interest  to  geophysicists  and  others  who  wish  to  determine  the  composition  of 
the  ocean  bottom.  The  plane  wave  reflection  coefficient  is  also  used  as  a  geometry  independent 
characterization  of  the  bottom.  As  such,  if  it  is  estimated  in  a  region  from  one  set  of  acoustic 
measurements,  then  the  fields  associated  with  an  arbitrary  source-receiver  geometry  in  that 
region  can  be  determined.  This  is  of  great  value  in  problems  of  acoustic  imaging. 

The  inversion  procedure  that  we  consider  in  this  chapter  was  proposed  by  Frisk, 
Oppenheim,  and  Martinez  [1  ].  It  requires  as  input,  coherent  measurements  of  the  pressure  field 
as  a  function  of  range  resulting  from  a  CW  point  source  over  a  horizontally  stratified  ocean  bot¬ 
tom.  From  this  the  (complex)  reflected  pressure  field,  Pg(r),  is  extracted.  The  Hankel 
transform  of  this  field  is  computed  to  provide  the  depth-dependent  Green's  function  as  a  func¬ 
tion  of  horizontal  wave  number:1 

Q(krf  so)  m  fr*  (r)J^kfr)rdr  (1) 

Finely,  the  plane  wave  reflection  coefficient  is  obtained  by  multiplying  the  Green's  function  by 

Ws  will  mbiHbu  ifcort—  Tdtrth-dspaadt  Or—n'i  function*  to  ’Om'i  htaettoa'. 


t 


terms  which  compensate  for  the  source  spectrum  and  for  the  source-receiver  distance: 

T(k, )  -  -i  Vki-kr2e~‘V^’i]‘  +,#IG  (*,  ,2  ,zq)  (2) 

This  entire  procedure  is  summarized  in  Figure  (V.  1.1). 

In  this  chaster  we  will  concentrate  on  the  estimation  oi  the  depth-dependent  Green’s  func¬ 
tion.  We  divide  the  issues  addressed  directly  into  the  categories  of  souroe-field  subtraction,  sam¬ 
pling,  windowing,  and  source- height  variation.  The  issue  of  source-field  subtraction  arises 
because  the  plane  wave  reflection  coefficient  is  directly  related  to  the  reflected  pressure  field  and 
not  the  total  pressure  field,  which  is  measured.  The  issue  of  sampling  covers  the  effects  caused  by 
having  die  pressure  field  available  for  computation  only  on  a  discrete  set  of  points.  We  discum 
both  the  effect  that  sampling  rate  has  on  the  estimate  of  the  depth-dependent  Green’s  function 
and  the  practical  problem  of  interpolation,  which  is  required  to  move  the  field  from  one  grid  to 
another  (often  to  compensate  for  missing  data  points).  We  develop  a  phase  unwrapping  pro¬ 
cedure  that  allows  us  to  interpolate  the  magnitude  and  unwrapped  phase,  which  vary  slowly 
compared  to  the  quadrature  components. 

In  the  section  on  windowing  we  discum  the  effect  that  having  the  pressure  field  available 
only  to  a  finite  range  has  on  die  estimate  of  the  depth  dependent  Green’s  function.  We  deter¬ 
mine  the  range  over  which  the  data  must  be  known  in  order  to  accurately  determine  die  depth- 
dependent  Green’s  function.  We  do  this  by  using  the  properties  of  the  Hankel  transform 
developed  in  Section  (IL6). 

In  the  section  on  source-height  variation  we  exploit  the  results  of  Section  (II. 6)  once  again, 
but  this  time  we  use  diem  to  determine  the  effect  that  variations  in  the  source-height  has  on  die 
estimate  far  die  depth-dependent  Green’s  function.  Such  variations  are  inevitable  during  the 
acquisition  of  real  data.  We  illustrate  these  effects  by  considering  the  effect  of  three  specific  van- 
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Figure  V.1.1  The  inversion  procedure  to  estimate  the  plane  wave  reflection  coefficient  from  the 
total  field  jenerated  by  a  CW  point  source  over  a  horizontally  stratified  bottom 
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V.2)  Sourct-fUld  subtraction 

la  this  section  we  consider  the  removal  of  the  source  field  when  the  source-receiver 
geometry  is  known.  In  Figure  (V.1.1)  we  showed  the  source  field  removed  before  die  Hank  el 
transform  We  did  this  because  conceptually  we  wish  to  deal  with  the  reflected  field  alone.  In  this 
section  we  show  that  numerically  it  is  better  to  remove  the  source  field  in  die  transform  domain, 
after  die  Hankel  transform  of  the  total  field  has  been  computed. 


Because  the  Hankel  transform  is  a  operator,  in  principal  the  estimate  for  the  Green's 
function  can  be  made  by  subtracting  the  source  field  either  before  transforming: 


G(kr,z,r0)  -  J 

o 


J^kfr)rdr 


or  by  subtracting  in  the  transform  domain: 


■  totVri+ii -st f 


G(kr,z,x0)  =  fPT(r)Ja(krr)rdr  -  - -J£k,r)rdr 

o  0  Vrz+(r-ro)z 

which  becomes  upon  performing  the  second  integral  analytically: 


(1) 


(2) 


G(k,jj o)  -  fP r(r)J0(krr)rdr  -  ^  k  -c‘  (3) 

If  Pj-(r)  is  available  only  over  die  finite  range  0<r  <rMt  then  the  field  integrals  can  only  be 
carried  out  to  rmn.  Substituting  rM  for  *  in  Equations  (1)  and  (3)  will  make  these  two  formu¬ 
lations  no  longer  equivalent  because  the  analytically  performed  integration  is  not  windowed. 

The  function  transformed  in  Equation  (1)  is  the  reflected  pressure  field,  Pg(r).  In  Section 
(IV. 3)  we  argued  that  the  reflected  field  decayed  asymptotically  as  y.1  If  the  total  field,  Pf(r), 

decays  asymptotically  faster  — ,  we  can  expect  that  the  formulation  of  Equation  (3)  to 

•aftttr  leas  from  windowing  effects  than  the  formulation  of  Equation  (1).  We  wiB  now  show  that 
die  total  field  does  in  fact  decay  faster  dun  the  reflected  field  alone.  In  fact,  by  transforming  the 

1)  la  tfes  abeam  at  topped  nodes. 
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total  pressure  field  and  then  subtracting  the  source  contribution  in  the  Hankel  domain,  we  have 

performed  the  dual  to  the  removal  of  the  — yJLmmm’  singularity  discussed  in  Section  (IV. 3a). 

V  k$  —kr 


Before  we  begin,  an  analogy  to  a  procedure  for  determining  die  Fourier  transform 

of  a  function  known  only  over  a  finite  range  but  with  a  large,  known  constant  offset  might  pro¬ 
vide  same  inright.  If  such  a  function  is  (Fourier)  transformed  directly,  die  offset  will  transform 
to  an  impulse  at  the  origin  which  is  smeared  into  die  rest  of  the  transform.  The  Mimaring  will 
occur  because  the  transform  is  taken  over  only  a  finite  aperture  (windowing).  The  alternative  is 
to  subtract  the  offset  from  the  function,  transform  the  result,  and  add  an  impulse  (with  a 
strength  which  is  determined  analytically  from  the  known  offset)  to  the  origin.  This  second  pro¬ 
cedure  will  give  superior  results  because  the  transform  of  the  offset  is  not  degraded. 


Transforming  the  reflected  field  alone  is  analogous  to  generating  die  Fourier  transform 
directly  from  the  the  field  frith  die  known  offset.  In  die  case  of  the  reflected  field,  however, 
instead  of  a  simple  constant  offset,  the  function  has  a  known  asymptotic  behavior.  It  decays  as 

y.  We  are  about  to  show  that  adding  the  source  field  to  the  reflected  pressure  field  is  analogous 
to  subtracting  the  offset  in  the  Fourier  transform  example.  In  the  Hankel  transform  case  we  are 
actually  considering  it  corresponds  to  subtracting  a  term  with  the  same  asymptotic  y  behavior. 

The  difference  will  decay  faster  than  — . 


We  begin  by  considering  the  Green’s  function  associated  with  the  total  field  for  xq>z 
which  is  given  by: 


v  —kr 


(4) 


The  ■— mw  term  is  the  source  term.  If  we  rewrite  Equation  (4)  in  terms  of  the  reflec¬ 


tion  coefficient  at  kg,  it  will  be  more  clear  why  adding  this  term  in  die  transform  domain 
corresponds  to  subtracting  the  asymptotic  behavior  in  die  pressure  domain.  We  must  use  dm  fact 
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that  for  all  bottom*  T(ho)  *  ~1  (except  the  degenerate  case  r(k,.)  ■  1  ).  We  write  Equation 
(4)  as: 

Gt(*„z,z0)  -  T7-T  ,  [r(i,)  -  r(*0)e‘aV*^O]e'V^^*°"l)  (5) 

Vktf-k2 

As  we  discussed  in  Section  (IV. 3)  the  asymptotic  behavior  in  the  pressure  domain  is  determined 
by  die  behavior  at  the  singularities  in  the  transform  domain.  [2  ]  At  kr  =  ko  the  phase  term, 

e  equals  1  so  that  unlike  the  Green’s  function  associate  with  the  reflected  field 

alone,  die  numerator  of  the  total  Green’s  function  approaches  aero  as  kr  approaches  the 

singularity  at  ko.  We  wish  now  to  determine  the  contribution  of  the  singularity  at 

kr  ■=  ko  in  die  total  Green’s  function  in  order  to  show  that  the  "softening"  introduced  by  the 
addition  of  the  source  term  has  made  the  associated  total  field  more  range  limited.  We  can 
bypass  a  great  number  of  issues  by  instead  considering  the  asymptotics  of  the  simplified  Green’s 
function: 

o)  -  - -5  [-1  ~  r(k0> ~*^ki }etVki (6) 

Vko  ~kr 

By  considering  Equation  (6)  we  exdude  those  issues  assodated  with  f(kr).  Our  examples 
in  the  synthetic  data  generation  section  showed  that  these  terms  do  not  give  rise  to  terms  which 

decay  as  slowly  as  y. 

Equation  (6)  is  the  Green’s  function  for  a  dipole  and  has  the  known  Hankel  transform: 

«o Vr* +(»-*,)* 

Pd(r)  *  *>  ,'  !  '  2  “  T7TT  »2  W 

Vr2+(z-r  o)2  Vr2+(x+z0)2 

It  is  well  known  that  tins  field  decays  asymptotically  as  -4-  and  that  this  asymptotic  behavior 

r* 

begins  more  quickly  when  zq  is  small  than  when  it  is  large. 

Since  die  total  field  will  be  more  range  limited  than  the  reflected  field,  it  is  better  to 
transform  the  total  field  numerically  and  subtract  the  (analytically  determined)  transform  of  the 


r 

i 

i 

c 
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incident  field  than  it  if  to  subtract  the  incident  field  before  transforming. 


V.3)  Sampling 
a)  Overview 

Typically,  data  is  not  available  on  the  grid  required  for  processing.  In  the  experiment  pro* 
Tiding  die  data  for  this  thesis,  for  example,  die  range  values  for  which  the  data  has  been 
obtained  is  determined  by  die  drift  rate  of  the  boat  and  the  source  away  from  die  receivers.  The 
individual  samples  do  not  occur  exactly  where  we  would  like  them  and  while  the  experiment  was 
designed  to  provide  samples  as  dose  to  the  Nyquist  rate  as  possible,  typically  the  number  of  sam¬ 
ples  on  averages  is  less  than  we  would  like.  Finally,  there  are  isolated  cases  of  missing  samples,  a 
reality  of  data  collection.  Section  (V.3b)  discusses  the  issues  associated  with  the  average  sam¬ 
pling  rate.  The  issues  associated  with  the  grid  in  general  are  discussed  in  Section  (V.3c). 


b)  Sampling  rate 

In  Chapter  (n.7)  we  saw  that  when  /  (r)  is  sampled  on  approximately  a  linear  grid  and  die 
transform: 


M 

F(p)  =  Sf(rV{hr)rdr  0) 

0 

is  computed  from  these  samples,  then  /  (r )  must  be  sampled  on  a  grid  at  least  as  fine  as  in 

order  to  correctly  perform  the  transform  for  F(p)  negligible  p>A .  In  this  chapter  we  consider 
the  transform  of  the  pressure  field,  to  obtain  the  depth -dependent  Green’s  function.  This 
transform  has  the  form: 


G(*„x+*0)  -  :7T'rlV^l,t,al  -  jfoC'VoCMWr  (2) 

V  0 

G(krj+z o)  is  negligible  for  kr>ko+*1  except  pomibly  near  die  poles  of  T(kr)  (for 


« 


1)  For  sows  tsiaBc>0. 


■nail  c)  because  when  kr  >kQ,  it  decays  exponentially.  Consequently  when  there  are  no  pedes  in 


IX*f)  for  real  kr,  then  the  pressure  field  need  only  be  available  on  a  grid  as  fine  as  with 

A 

A  *  i0+e  to  accurately  determine  G  (Jkr  ,z  +zq)  in  the  region  0 <kr  <ko .  If  we  wish  to  obtain 
G(kr,z  +*o)  *»  region  where  it  is  exponentially  decreasing  (kr>kd),  however,  we  must  sam¬ 
ple  fast  enough  to  represent  die  signal  in  that  region  as  well. 

If  a  pale  is  present  in  T(kr)  at  kr  =  k,t  the  Green's  function  will  be  significant  near  kr( 

despite  the  exponential  decay.  If  the  presence  of  the  pole  is  ignored  and  the  field  is  transformed 
it  ir 

on  the  grid  -j-~,  then  the  pole  will  be  aliased  into  the  Green’s  function  at  lower  wave  numbers. 

If  there  is  only  one  pole  present  we  can  write  die  Green’s  function,  G  (kr)  (we  suppress  die 
z  variation)  as: 


G(*,)  «<?(*,)  + 


o-i 


for  kr>k 


(3) 


K 

*  j  * 


The  results  of  Chapter  (II.7)  show  that  if  the  pressure  field  is  transformed  on  the  grid  — ,  then 


the  aliased  Green’s  function  computed  will  approximately  be  given  in  the  region  0<ir<Jt  by: 


C(kr)  =  G(kr)  -  V2k/k,-lG(2k  ~kr)  =  G(kr)  ~  V2k/kr-l  -  (4) 

(2k-kr)2-krJl 

so  that  the  Greens  function  at  2k^—kr.  wifi  be  corrupted. 

If  the  amplitude,  o_ j,  is  very  (which  would  be  the  case  for  large  source-receiver 
geometries)  and  some  smearing  is  present  due  to  windowing  (the  field  is  not  measured  out  to 
ranges  where  the  trapped  mode  dominates),  we  may  not  see  the  pole’s  effect  and  it  can  safely  be 
ignored. 

In  general  the  possibility  of  trapped  modes  must  be  considered  before  deciding  upon  a  sam¬ 
pling  rate,  particularly  in  geometries  with  small  source-receiver  heights.  For  such  geometries  it  is 

not  always  sufficient  to  sample  at 


c)  Sampling  grid 


When  data  is  not  available  on  the  grid  required  for  processing  we  must  first  interpolate. 
Successful  interpolation  is  possible  only  if  the  signal  is  adequately  represented  by  the  original 
samples.  If  we  know  only  that  our  signal  has  a  Hankel  transform  which  is  negligible  beyond 


some  bandwidth,  A,  then  the  signal  is  adequately  represented  by  samples  on  the  grid  —  for 


C  iA  and  where  X,  n*0,l,2,  *  *  *  are  the  ordered  zeros  of  [2  ]  This  is  true  in 

theory.  In  practice,  if  the  the  samples  are  not  originally  spaced  as  required,  it  may  be  impossible 
to  actually  perform  the  interpolation  onto  another  grid.  If  the  samples  are  only  available  on  the 


grid  —  with  C  <A ,  then  it  is  not  possible  to  interpolate  without  making  additional  assump- 

V* 


tions. 


We  will  show  that  for  the  das  of  pressure  fields  examined,  an  additional  assumption  seems 
reasonable.  This  assumption  makes  interpolation  posable  even  when  the  sampling  rate  is  slightly 
too  low.  We  will  assume  that  the  magnitude  and  phase  of  our  pressure  fields  are  smooth  com¬ 
pared  to  their  quadrature  components.  Figure  (V.3c.l)  shows  the  magnitude  of  the  pressure  field 
associated  with  a  point  source  over  a  pressure  release  bottom.  Between  calculated  points  the 
curve  varies  so  little  and  so  regularly  that  a  plot  of  the  points  appears  to  be  a  smooth  line.  Fig¬ 
ure  (V.3c.2)  shows  the  result  of  first  subsampling  the  points  plotted  in  Figure  (V.3c.l)  (which 

K 

were  available  on  the  grid  )  by  a  factor  of  two,  and  then  interpolating  back  onto  the  original 

grid  using  splines.1  The  differences  between  die  two  curves  are  negligible. 

We  can  compare  this  succeariul  interpolation  to  the  result  of  subsampling  the  quadrature 
components,  spline  interpolating,  and  computing  the  magnitude.  The  result  of  this  operation  is 
shown  in  Figure  (V.3c.3).  The  apparently  smooth  line  comes  from  die  subsampled  set  of  values 

which  die  splines  was  constrained  to  match  in  the  quadrature  components.  It  actually  consists  of 

U  Spfesa  wan  used  beauts  We  original  grid  Is  gtvsa  by  X.,M  ,  wbsrt  X,  art  the  order  ad  tarot  of 
J o(x  ).  TUs  grid  is  usovaa  ad  makst  otbsr  interpolation  schema  last  desirable. 


MAGNITUDE 


•“  A 


A 


Figure  V.3c.2  The  result  of  sub  samp  ling  and  interpolating  the  magnitude  function  sho 
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every  other  point  of  the  displayed  curve.  The  surrounding  scatter  is  the  magnitude  of  interpo¬ 
lated  points  supplied  by  the  splines.  Clearly,  splines  did  not  successfully  interpolate  the  quadra¬ 
ture  components. 

Figure  (V.3c.4)  shows  the  phase  of  Figure  (V.3c.l)  computed  by 

4  -  tan'1  ■£-  (1) 

rT 

where  Pt  is  the  imaginary  component  of  the  field  and  Pr  is  the  real  component.  The  rapid  vari¬ 
ation  in  0  suggests  that  it  is  not  adequately  represented  by  the  grid  upon  which  it  is  presented.  § 
is  not  the  only  representation  of  the  phase  of  the  pressure  field,  however. 

d)  Unwrapping  the  Phase 

The  phase  displayed  in  Figure  (V.3c.4)  is  the  principal  value  of  the  phase,  often  referred  to 
as  the  wrapped  phase.  The  wrapping  comes  about  because  of  die  ambiguity  concerning  which 
phase  should  be  associated  with  the  quadrature  components.  If  0  satisfies: 

Me1*  «  f»,  +  i  Pi  (1) 

then  so  must  0  •+•  2-irm  where  m  is  any  integer,  since 

Afe'f*  +  2vl  •>  =  Me‘*e2v‘  m  =  Me1*  =  Pr  +  iPt  (2) 

Given  just  Pr  and  Pi  there  is  no  way  to  determine  the  correct  value  of  m .  The  arctan  routine 

used  by  Fortran  follows  the  convention  of  choosing  m  such  that 

— ir  <  §  ■  0  +  2‘irm  <  ir  (3) 

The  output  value  0  is  the  principal  value  of  the  phase,  or  the  wrapped  phase. 

If  the  phase  of  the  pressure  field  were  approximately  increasing  at  a  rate  of  k(/t  where 

n-V^TcT  —  zq)2  and  the  field  were  sampled  at  the  Nyquist  rate  of  1/2*0,  then  the  phase 

difference  between  samples  would  be  roughly  it  and  the  wrapped  phase  every  sample  or  two 

would  suffer  a  jump  to  a  different  m  in  order  to  satisfy  the  condition  -it  <  0+2irm  <  sr. 

This  would  obscure  any  underlying  regular  behavior  that  we  expect  from  most  phyacal 

phenomena.  Tbeee  frequent  jumps  are  responsible  for  the  rapid  oscillation  apparent  in  the  phase 


of  Figure  (V.3c.4).  We  wish  to  compensate  for  the  wrapping  that  takes  place  when  the  principal 
value  of  the  phase  is  generated.  To  do  so  we  must  make  use  of  our  knowledge  of  how  the  phase 
is  varying  from  point  to  point. 

We  conjecture  that  the  phase  of  the  total  field  is  dominated  by  a  component  at  the  water 
wave  number  associated  with  the  dominant  specular  path  and  that  the  remaining  portion  of  die 
phase  is  slowly  varying  compared  to  the  —mpltwg  rate.  We  write 

PT(r)-M(r)e‘ (4) 

where  M  (r)  and  0(r)  are  real,  and  write 

«(r)  =  M+e(r)  (5) 

where  R  -  Vr2  +  We  will  call  k<)R  the  modeled  portion  of  the  phase  and  e(r)  the 

residual  phase.  We  are  going  to  show  that  as  long  as  the  residual  phase  is  sampled  fast  enough, 

we  can  reconstruct  the  true  phase. 

In  this  notation  the  difference  in  true  phase  value  from  sample  to  sample  can  be  written: 

»(%)  "  <K r.-x)  -  k0(R„  -  +  €(*„)  -  *(*„_!)  (6) 

so  that 

0(0  -  e(rM-J  -  k0(R„  -  *,-!>  *  e(R.)  -  «(*„-!)  (7) 

Precisely  stated,  our  requirement  that  the  residual  phase  be  slowly  varying  compared  to  the  sam¬ 
pling  rate  is: 

l«(*J  -  «(K«-l)l  <7r  for  all  R„  (8) 

To  unwrap  the  phase  we  first  form: 

6(r.)  "  "  k0(R,  -  H.-O 

-  9(r„)  -  m„2ir  -  8(%-i)  -  m.-i2ir  -  k^R,  -  R„-i) 
m  9(r„)  ~  e(r,_!)  -  2ir(ma  -  ma.x) 

■  «(*.)  ” 

from  the  measured  data.  We  now  do  the  unwrapping  by  <***■*■!  m0  ■  0  and  picking  the 
integers,  m„,  (it  ■  1,  2,  •  •  •  )  sequentially  to  satisfy: 


1 6(0  -  -  k0(Rm  -  *„_!>  -  2ir(m(l  -  m.-O!  =  |«(*.)  -  «( Rm.{)\ <it  (10) 

and  define  the  unwrapped  phaie  to  be: 

8(%)  *  Hr*)  +  (11) 

Figure  (V.3d.l)  shows  the  result  of  running  this  algorithm  on  the  phase  of  the  synthetic 
data  with  magnitude  shown  in  Figure  (V.3c.l)  and  wrapped  phase  shown  in  Figure  (V.3c.4). 
The  resulting  phase  is  dominated  by  die  linear  term  kg R  we  defined  in  our  model.  Figure 
(V.3d.2)  shows  the  residual  phase.  The  smooth  and  small  variation  of  the  residual  phase  over 
die  intervals  ,  rm  j  for  all  n ,  is  a  strong  confirmation  of  our  phase  unwrapping  assump¬ 

tion. 

Figures  (V.3d.3a)  and  (V.3d.3b)  present  the  magnitude  and  residual  phase  of  the  fast  bot¬ 
tom  example  of  Section  (IV. 3a).  For  this  example,  too,  the  residual  phase  is  well  behaved. 

Figures  (V.3d.4a)  and  (V.3d.4b)  present  the  magnitude  and  residual  phase  of  the  slow 
speed  layer  example  of  Section  (IV.3b).  For  this  example,  too,  the  residual  phase  is  well 
behaved.  The  field  in  this  example  was  shown  to  be  dominated  in  the  far  field  by  the  contribu¬ 
tion  due  to  the  pole  beyond  the  water  wave  number.  The  upward  slope  of  the  residual  phase 
apparent  in  Figure  (V.3d.4b)  reflects  the  fact  that  this  pole  is  contributing  the  dominant  com¬ 
ponent  to  the  phase  (in  the  far  field)  which  is  slighdy  larger  that  the  kgR  term  subtracted  out. 


e)  Interpolating  the  magnitude  and  unwrapped  phase 

In  Figures  (V.3c.l)  and  (V.3c.2)  we  showed  that  the  magnitude  of  the  dipole  field  could  be 
up-sampled  from  the  grid  —  to  In  Section  (V.3d)  we  saw  that  the  unwrapped  and  resi- 
dual  phase  components  enjoy  smooth,  regular  variation  ideally  suited  for  interpolation.  Figure 
(V.3e.l)  shows  tile  residual  phase  for  the  dipole  field  up- sampled  from  the  —  to  the  grid. 

X  me 


It  is  indistinguishable  from  the  residual  phase  generated  on  tire  ~  grid  shown  in  Figure 
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Figure  V.3d.l  Unwrapped  phase  corresponding  to  Figure  V.3c.4 
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Figure  >'.3d.2  Residual  phase  corresponding  to  Figure  V .3c.4 
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Figure  V.3d.4b  Residual  phase  for  slow  bottom  example. 


PHASE 


We  now  show  that  for  the  dipole  field  we  can  actually  interpolate  the  magnitude  and 
unwrapped  phase  to  increase  the  sampling  rate  of  the  quadrature  components  of  the  field.  This 
allows  us  to  determine  the  Green’s  function  to  a  higher  horizontal  wave  number  than  the 
Nyquist  criteria  applied  to  our  original  sampling  scheme  would  have  us  believe. 


We  recall  from  Section  (Q.7)  that  if  the  pressure  field  for  the  dipole  on  the  grid  were 

transformed  and  displayed  in  die  range  0<Jk,<2  the  result  would  be  severely  aliased  and  com¬ 
pletely  inaccurate  in  the  region  l<tf  <2.  To  obtain  a  transform  accurate  on  0<ir<2,  the  qua¬ 
il,, 

drature  components  must  be  at  least  sampled  on  the  grid  We  can  still  obtain  the  transform 


in  the  range  0<kr<2,  never-the-less,  by  interpolating  the  field  onto  the  grid  through  its 


magnitude  and  unwrapped  phase.  Figures  (V.3e.2a)  and  (V.3e.2b)  show  the  magnitude  and 
phase  of  the  transform  generated  by  such  a  procedure.  First  the  magnitude  and  residual  phase  of 


p(-y-)  were  generated.  These  were  up-sampled  as  just  discussed.  From  this  up-sampled  magni¬ 


tude  and  residual  phase  (and  the  modeled,  kgR  ,  portion  of  the  phase)  the  associated  quadrature 
components  were  generated.  This  was  transformed.  Figures  (V.3e.3a)  and  (V.3e.3b)  show  the 


magnitude  and  phase  of  the  Hank  cl  transform  of  p  (-j— )  generated  without  interpolation.  Only 


small  differences  in  the  magnitude  are  apparent.  The  phase  curves  also  display  only  small  differ¬ 
ences  though  in  the  inhomogeneous  region  (where  the  phase  is  oscillating  rapidly  as  evidenced  by 
the  two  parallel  lines)  die  small  difference  has  caused  a  slightly  different  picture  of  the  oscilla¬ 
tions.  By  contrast.  Figures  (V.3e.4a)  and  (V.3e.4b)  present  the  magnitude  and  phase  of  the 


Hankel  transform  of  p(— )  up-sampled  by  direct  spline  interpolation  of  its  quadrature  com¬ 


ponents.  Clearly,  once  again,  a  direct  interpolation  of  the  quadrature  components  did  not  work. 

We  apply  this  scheme  for  interpolating  the  magnitude  and  residual  phase  to  the  field  of  the 
fast  bottom  example  of  Section  (TV. 3a).  We  first  generate  the  magnitude  and  unwrapped  phase 


PHASE 


Flgnre  V.3e.2b  Phase  of  Han  ltd  transform  of  up-sampled  field  for  pressure  release  example 
shown  in  Figures  V.3e.2  and  V.3e.l 
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Figure  V.3*.3a  Log-magnitude  of  Hankel  transform  of  field  for  pressure  release  example 
without  up-sampling 
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Figure  V.3e.3b  Phase  of  Hank  el  transform  of  Held  for  pressure  release  example  without  up- 
sampling 
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Figure  VJi.4b  Rum  of  Hiakel  transform  of  field  for  the  pressure  release  example  that  has 
been  op-sampled  through  its  quadrature  eompoomits 
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of  1(2^3*)’  interpolate  up  to  the  grid  generate  the  quadrature  components,  and  then 
Hankel  transform.  The  magnitude  and  phase  of  the  result  is  shown  in  Figures  (V.3e.5a)  and 
(V.3e.5b).  The  magnitude  and  phase  of  the  correct  transform  (of  generated  without 

/a 


using  this  interpolation  scheme)  is  shown  in  Figures  (V.3e.6a)  and  (V.3e.6b).  We  see  that  the 
Hankel  transform  of  the  up-mmpled  data  and  the  Hankel  transform  of  the  data  originally  avail¬ 
able  on  the  fine  grid  do  not  agree  exactly.  Figures  (V.3e.7a)  and  (V.3e.7b)  present  the  magni¬ 
tude  and  phase  of  their  complex  difference  and  Figure  (V.3e.8)  presents  the  magnitude  of  the 
Hankel  transform  of  that  complex  difference.  This  transform  represents  the  errors  made  in  the 
pressure  domain  by  our  up-sampling  procedure  that  gave  rise  to  the  error  in  the  Green's  func¬ 
tion.  We  see  that  practically  all  the  error  energy  was  concentrated  at  the  origin.  This  error  could 
be  due  to  a  breakdown  in  our  phase  unwrapping  aammption  near  the  origin  or  to  a  poor  han¬ 
dling  of  the  rapid  change  in  magnitude  by  the  splines.  This  problem  can  be  correctedby  a  dense 
sampling  of  the  original  field  near  the  origin  so  that  there  is  no  room  for  interpolation  error 


thqpe. 


f)  Phase  unwrapping  errors 


At  this  point  we  consider  briefly  the  kinds  of  error  that  might  be  expected  when  die 
assumption  underlying  this  phase  unwrapping  technique  is.  violated.  If  for  some  n 


!»(*,,)  -  efo.x)  -  to (Re  -  Re- 1) I  >  *  (D 

the  wrong  ma  will  be  chosen.  From  that  point  on,  each  mt  (!■«,«  4-1,  •••  )  chosen  by  the 

procedure  will  also  be  wrong  by  die  same  amount  A  plot  of  this  error  is  a  step  function  of 

height  -  mm  centered  at  n  as  shown  in  Figure  (V.3f.l).  If  multiple  violations  occur,  the 

error  will  look  like  the  sum  of  nsp  functions  as  illustrated  in  Figure  (V.3L2).  The  smoothness 

apparent  in  the  residual  phase  in  all  of  our  examples  suggests  that  no  errors  have  occurred. 

B  the  phase  unwrapping  scheme  is  used  to  interpolate  dm  field,  these  errors  are  not  serious. 


As  part  of  the  interpolation  procedure,  the  quadrature  components  are  regenerated  from  the 
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Figure  V  J*.5»  Log-magnitude  of  Hinkel  transform  of  the  field  for  the  fast  bottom  example 
after  it  has  been  up-sampled  through  its  magnitude  and  residual  phase 
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Figure  V.3*.6a  Log-magnitndc  of  Hankel  trmcsfonn  of  the  field  lot  the  fut  bottom  example 
generated  originally  on  a  fine  grid 
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Figure  V.3t.7b  Phase  of  the  complex  error  in  the  Hankel  transform  of  the  up-sampled  field  for 
die  fast  bottom  example 
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interpolated  magnitude  and  phase.  The  error  curve  shown  in  Figure  (V.3f.2)  would  have  no 
affect  on  the  quadrature  components  regenerated  from  the  unwrapped  phase.  In  general,  after 
the  interpolation,  the  error  curve  will  not  be  a  simple  sum  of  steps  but  will  be  smeared  by  the 
interpolator.  This  will  usually  affect  the  quadrature  components.  If  the  interpolator  is  well 
chosen,  the  leakage  will  be  small  and  ««»»»«*»  to  the  area  near  die  error.  Finally  we  note  that 
errors  in  the  phase  unwrapping  scheme  will  occur  when  the  unmodeled  portion  of  the  phase  is 
varying  rapidly  between  samples.  When  this  happens  the  interpolator  is  likely  to  have  difficulties 
even  without  errors  in  the  unwrapped  phase  and  this  scheme  is  probably  not  appropriate. 

V.4)  Windowing 

In  Section  (II.6)  we  stated  that  in  terms  of  resolution  the  Hankel  transform  behaves  very 
much  like  a  Fourier  transform.  We  wish  to  consider  the  resolution  required  to  generate  the 
Green's  function  and  the  window  that  this  implies. 

The  total  Green's  function  is  given  by 

Gr(*,,r)  *  y  when  r>,0  (1) 

vk$  ~kr 

The  most  rapid  variations  in  Gj  (excluding  possible  poles  in  the  reflection  coefficient  beyond  the 

i^/ti -*r*(f+**> 

water  wave  number)  occur  near  kr  “  k0.  When  T(kr)  is  smooth  compared  to 

-*,*(*  +*#) 

the  rapidity  at  these  variations  is  dominated  by  the  g  y  <■  term.  With  a  windowed 

sample  of  the  pressure  field  we  can  not  hope  to  determine  die  exact  behavior  of  G  at  kr  m  Iq, 

when  z  +z0  is  large  the  rapid  variation  in  G  is  due  to  primarily  the  term.  We 

can  obtain  an  ad  hoc  estimate  of  the  resolution  we  require  by  considering  the  lob*  widths  associ¬ 
ated  with  die  phase  for  kr  near  Iq* 


That  is  we  define  kr  ^  by  the  relation: 


[ela] 


and  define  the  lobe,  8a  by: 


V*0  +*o)  ■ 


:  we  use  (*<>“*#•,«  )(*o+*r,«)  *  (  ■)2  then  when  kftlt  ~  k  we  have: 


=  57  7% 


s-  -  37  '  37 

Section  (Q.6c)  indicates  that  the  required  window  width,  B ,  is  related  to  the  desired  resolution 


approximately  as: 


w  .3  '  J_  _  3(x-t»z0)2  i _ 1_ 


&  8. 


w2  2i0  2n  — 1 


Thus  to  resolve  the  lobe  doses  to  k  when  z  +z<)  *  136  and  *o  *  .9246  we  require  a  win¬ 


dow  of  about: 


-3,(1013>-  3040 


meters 


V.5)  Source-Height  Variation 
a)  General  expression 

The  procedure  proposed  to  estimate  die  plane  wave  reflection  coefficient,  r(tr),  and 
shown  in  Figure  (V.1.1),  requires  that  die  pressure  field  be  measured  with  the  source  at  a  fixed 
height,  zo-  [1  ]  Frequently,  experimental  conditions  cause  the  source-height  to  vary.  In  this  sec¬ 
tion  we  will  explore  the  effect  that  a  varying  source-height  has  on  die  estimate,  f  (*f). 

q|  considering  the  effect  of  a  varying  source-height  on  the  estimate  for  the  plane 
wave  reflection  coeffident  directly,  we  will  consider  its  effect  on  the  depth- dependent  Green’s 


function  given  by:1 


X)  Wawffl 


Oram’s  function  cm  Z  sad  Xq. 
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G(M  -  fr,(r?y^k,,)rdr  -  «> 

0  V  Kq  ™ 


(i) 


which  is  die  Hankel  transform  of  the  reflected  pressure  field.  By  considering  the  effect  on  the 
depth-dependent  Green’s  function  we  can  make  use  of  the  properties  of  the  Hankel  transform 
that  we  derived  in  Chapter  n.  The  plane  wave  reflection  coefficient  is  determined  by  multiplying 
the  estimate  of  the  depth-dependent  Green’s  function  by  terms  which  compensate  for  the  source 
strength  and  the  source-receiver  separation  as  was  shown  in  Equation  (V.1.2). 

We  consider  the  effect  of  a  source  height  given  by 


*0(0  “  *o+*(0  (2) 

To  explore  the  effect  of  Hankel  transforming  a  pressure  field  measured  at  a  source  height  that  is 

a  function  of  range,  we  write  the  Green’s  function  estimated  by  Hankel  transforming  this  field 


as: 


m 

m£pK{rj(r))J0(k,r)rdr 

We  now  define: 


9 

m  iJeiVki~^r)J^)J0(krr)rdr  (4) 

0 

which  with  (3)  becomes: 

9 

6(.k,)  -  fawfoMt  (5) 

Equations  (4)  and  (5)  exactly  describe  the  effect  that  source-height  variation  has  on  the 
estimate  of  the  depth-dependent  Green's  function.  As  written,  however,  they  do  not  provide 
much  insight  into  what  variations  are  tolerable  or  into  the  qualitative  effect  of  source-height  vari¬ 
ation.  To  provide  us  with  this  insight  we  develop  an  approximate  expression  for  H  (kr  ,£)  by 
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using  the  windowing  result  of  Section  (D.6).  We  do  this  by  considering  &h(r)  u  a  wja. 

dow.  The  result  will  be  reasonable  provided  that  the  Fourier  transform  (in  r)  of  &*(')  ^ 

narrow,  as  discussed  in  Section  (II. 6b). 

We  write 


Bfa,Q  -  tfel'*(r)Jofr\J0(krr)rdr  with  ««  \402-*2 
o 

,  v 

The  Hankel  transform  of  Jq(v)  equals  — - - so  that : 


(«) 


(Vfc7-&fe°)*ws(  kr) 


where 


(7) 


W{(*r )  ■  / (8) 

— ao 

This  provides  us  with  approximate  expressions  for  the  kernel,  H(kr ,£),  and  an  approximate 
expression  for  the  estimated  Green’s  function  in  terms  of  the  actual  Green’s  function: 


H  fa  .€)  *  ~^j^W  (fa  ”0 

M 

6  fa)  * 

In  the  following  sections  we  apply  this  result  to  some  special  cases. 


(9) 


b)  Particular  variation 
I)  No  xourct-htight  variation 

When  the  source-height  is  constant,  h(r)  ■  0.  For  this  case  our  approximate  result  above 
gives  Wfa-Q  -  2w8(*r-0  and  <$(*,)  *  G  (kr ),  which  is  u  we  would  expect 

* 

D  Linear  tourct-htight  variation 

If  the  source-height  varies  linearly  then  h(r)  m  ar  and 
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Wk(kr)  =  Jeilua+k,]rdr  =  2Tth(k,+aVk$-?)\  This  givet  ut 
— » 

X  _______ 

G{kr)  -  ^5^°  (*)»(*, -* +«  Vkl-?)di  (1) 

To  evaluate  the  integral  we  have  to  amplify  the  argument  of  the  delta  function.  We  define 
«  ■  so  that 

c  _  a  +  Vs2-(l+a2)(s2-o2k§) 

*  l+o2  1 ' 

Substituting  into  Equation  (1)  we  have: 


tv,)  =  v-Z^swo 


f  —ds 


2o2+£-j 


Where 


*r  +  Vi  ^-(l+a^^-o^o  ) 

*° - l+o  2 

*,+oV(o2+l)*|-*,2 

o2+l 


assuming  that  fc,  »  real.  We  see  that  <$(kr)  is  a  distorted  verson  of 
C  |^  +  - — ^ L |+^-° — —  .  This  approximate  analysis  also  correctly  indicates  that  at  die 

l  a 2+l 

slope  of  die  linear  variation,  a ,  goes  to  zero  6  (kr)  goes  to  G  ( kr ). 


Hi)  Sinusoidal  source-height  variation 

When  the  source-height  variation  is  given  by 

*(') 


. )  .  .  j,, 


w,(* 


1)  Ptoridad  {  <  Jk .  Tbs  iasgral  It  aot  dafiaad  ter  goapisx  arfuanta  sad  a  dUhnat  aalyds  woald  b# 
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with  <i>  *  a  -£2- 


To  perform  the  integral  in  Equation  (2)  we  eipand  the  exponential  in  its  Taylor  aeries  to 
obtain: 


-0  "  ■ 
%  m  I  J 

1*0  n  1  — 3T 


(3) 


2ir  £  -^^8(*r-«o) 

■  »o  *  • 


We  use  this  to  determine  the  effect  of  the  real  cosinusoidal  variation: 

Mlv  x  nM  -lor 

Hr)  -  - . -  (4) 

Substituting  (4)  into  Equation  (V.5a.8),  for  cosinusoidal  variation  W{(tr)  is  seen  to  be  given  by: 

W€(*r)  -  fe  2  e  2  e*'rdr  (5) 

-3B 

Equation  (S)  is  the  Fourier  transform  of  a  product  of  terms  in  the  form  of  Equation  (2).  Conse¬ 
quently  we  can  write  Equation  (5)  as  the  convolution  of  terms  in  the  form  of  Equation  (3): 


W&r) 


2ir 


2*2  ii^J-b(kr-na) 

mm0  *  ' 


2*2  2 

a  -0»  -0 


iLs)L 


* 


2TT2&£-6(kr+na) 
a  “0  **' 


(6) 


n!m ! 


S(Jtr  ~(n  -m)a) 


II  we  perform  the  integration  (V.Sa.9)  we  obtain: 


<?(*,)* 


1 


*f  >(»  -m)« 


•  tail 


-m)aC  (*,-(«  -w)«) 


(7) 


The  cosiaiisoidal  source-height  variation  with  an  amplitude,  a ,  and  a  frequency  a,  has  the 
effect  of  reverberating,  or  comb  filtering,  y/k^G(kr)  in  two  dimensions.  The  impulses  of  die 

raqolrad  tor  tUi 


filter  are  (paced  at  the  variation  frequency  a.  The  weightings,  which  we  define  to  be 
w  (kryn  ,h),  determine  the  envelope  of  the  reverberation.  They  depend  on  die  amplitude,  a : 


w(kr\mjt) 


We  can  write  the  estimate  for  die  Green’s  function  in  terms  of  these  weighting  functions  as: 


i  22 

<?(*,)«  -^r  kr>(H-m)a  " (*, -(*  -m )*j* ,n)Vk,  -(*  -m )a G (*, -(»  -m )a )  (9) 

The  weighting  functions,  w (kr  -  (n  —m) a‘/n  ,n),  are  greatest  when  m  «  n  *  nma  and 
decay  rapidly  from  that  point  in  m  and  it .  This  result  can  be  shown  by  replacing  the  factorials  in 

(8)  with  Stirling’s  approximation  (excellent  even  for  small  n :  n !  ~  jyj  )  and  defining 


la  Vkj  ~k? 


.  The  weighting  functions  then  become: 


„(*,*.  ^)-44= ^4-kr  -4-f^v  do, 

v  ’  ml  n\  V2tr m  ( m  J  ^2 irn  [  a  J 

(—  [  term  has  its  maximum  at  m  «x  and  falls  off  in  m  with  greater  than  geometric 
m  ) 


decay.  The  .  r—  term  pulls  this  maximum  only  very  slightly  lower. 

vm 


The  remit  is  that  w(kryn  ,*)  is  large  for  m ,  n  ~~ ‘S/kfi-k?  and  mall  elsewhere.  When 
/ kft  —k?a  S,  ® -p  J  «1  we  can  ignore  die  (n  -m  )a  term  in  (9)  so  that  (j  (kr)  is  given' 


approximately  by 


.  1  SI 

<2(*r)  83  kr  >(a  -si 


W-k? 


fk,-(n  -m)aG  (*,-(»-»»)«) 


By  defining: 
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c(*,*)-(i)*  2 

mmmax{ 0,-») 


(n+m)\ 


we  can  obtain  the  result,  valid  for  «  1  that 


2  C (k, )Vkr-naG (kr -n a)  (13) 

ft  -  -» 

A  perspective  plot  of  C(krtn)/(i)*  is  presented  in  Figure  (V.Sc.la)  for  the  case  a  =  3 
and  &0  ~  -9246159.  The  back  of  this  figure  corresponds  to  kr  =  Aq  and  consequently  1,  *  0. 
The  Green's  function  in  this  region  corresponds  to  plane  wave  components  of  die  field  which  are 
directed  entirely  in  die  radial  direction  and  which  do  not  vary  in  z.  Figure  (V.5c.lb)  presents 
the  slice  of  Figure  (V.5c.la)  corresponding  to  this  region,  C  (k0,n).  C  (k0,n)  is  zero  everywhere 
except  at  n  =  0,  where  it  is  1.  Referring  to  Equation  (13)  we  see  that  the  degraded  estimate  of 
the  Green’s  function  at  kr  *  k0  is  given  by: 


£(±o)  *  -T7T  2  Cik^Wko-ncxGik.-na)  (14) 

Substituting  for  C  (ko,n  )  in  Equation  (12)  we  see  that 

<?(*o)-G(*0)  (15) 

The  portion  of  the  spectrum,  kr  ~ko,  corresponds  to  field  components  that  do  not  vary  in 

x .  It  Is  reasonable,  then,  that  the  cotinusoidal  source-height  variation  did  not  affect  that  portion 
of  the  angular  spectrum. 

hi  Figure  (V.Sc.la)  moving  forward  towards  the  leading  edge  corresponds  to  decreasing  kr 
and  increasing  k, .  With  deer  eating  kr,  C(krTn)  becomes  increatingly  less  impulsive,  indicating 
greater  amounts  of  degradation.  Figure  (V.5c.lc)  presents  the  slice  C(0^i).  This  time 


corresponds  to  that  portion  of  foe  angular  yectrum  which  has  the  maximum  amount  of  vertical 
variation.  In  Figure  (V.Sc.lc)  the  value  C(0,0)  is  not  even  as  large  as  foe  adjacent  values, 


O 
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Figure  V.Sc.ia  Perspective  plot  of  C  (fcr  ,n )/(i )"  with  a  =3 


C(krtn) 
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C  (0,-1)  and  C  (0,1).  The  Greens  function  will  be  degraded  by  cosinusoidal  source-height  vari¬ 
ation  in  this  region. 

Figure  (V.Sc.2)  presents  a  perspective  plot  of  C(Jfcr,n)/(i)*  for  the  case  a  =  12  and 
k0  =  .9246159.  Once  again  C(k0,n)  is  the  discrete  delta  function,  8(/t),  and  the  Green's  func¬ 
tion  wiD  not  be  degraded  at  kr  »  k(y  Because  a  is  larger  now,  ■— VjfcJ-fc]?  of  Equation  (12) 

grows  more  rapidly  as  kr  becomes  smaller  than  it  did  for  a  *  3.  As  a  result  the  figure  shows 
that  serious  degradation  begins  fox  kr  much  closer  to  Jko-  The  increased  amplitude,  a,  has 
resulted  in  an  increased  amount  of  degradation.  The  product,  ak,  =  a  Wk^—k2',  determines 
the  severity  of  this  effect. 

We  note  also  that  because  of  the  (i)"  factor  in  Equation  (12),  the  phase  of  C (Jtr,n) 


ically  affect  the  phase  of  the  estimated  Green's  function,  d(kr),  even  before  it  significantly 
affects  the  magnitude. 

Thus  we  have  seen  that  the  effect  of  sinusoidal  source-height  variation  is  to  comb-filter  the 
estimate  of  y/k^G(kr).  The  spacing  between  impulses  in  the  comb  filter  is  the  frequency  of  the 
source-height  variation.  The  amplitude  of  the  source-height  variation  and  the  vertical  wave 
number,  to  ~kr2,  together  determine  the  weightings  of  the  impulses.  When  the  product  of  the 
amplitude  and  the  vertical  wave  number  is  small,  the  only  contribution  comes  from  the  low  lag 
components.  As  this  product  increases,  the  higher  lag  components  begin  to  contribute  and  the 
comb-  filtering  will  become  increasingly  apparent.  If  the  frequency  of  the  source-height  variation 
is  very  low,  causing  the  spacing  of  the  impulses  in  the  filter  to  be  very  small,  the  degradation 
will  appear  as  a  smearing. 

V.S)  Summary 

In  this  chapter  we  have  studied  the  iames  associated  with  the  inversion  of  pressure  field 


L 
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data  through  the  Hankel  transform  to  estimate  the  depth-dependent  Green's  function  and  the 
plane  wave  reflection  coefficient.  We  have  developed  a  phase  unwrapping  procedure  that  allows 
us  to  interpolate  the  magnitude  and  unwrapped  phase  and  thereby  determine  from  the  set  of 
field  samples  available,  the  values  of  the  field  at  the  ranges  we  require  for  processing.  We  have 
also  shown  that  it  is  better  to  estimate  the  total  depth-dependent  Green’s  function  from  the 
Hankel  transform  of  the  total  field,  and  to  later  remove  die  affects  of  the  source.  Finally,  we 
have  examined  the  effects  of  source-height  variation  to  help  us  understand  the  possible  degrada¬ 
tion  that  this  effect  would  would  introduce  into  the  depth-dependent  Green’s  function  estimated 
from  real  data. 

We  are  ready  to  perform  a  preliminary  processing  of  real  data. 
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CHAPTER  VS) 

Inverting  Real  and  Realistic  Data 

In  Chapter  V  we  described  the  procedure  for  inverting  coherent  field  measurements  arising 
in  response  to  a  CW  point  source,  to  obtain  the  plane  wave  reflection  coefficient.  In  that 
chapter  we  addressed  some  of  the  practical  issues  that  must  be  faced  when  real  data  is  to  be 
inverted.  In  this  chapter  we  perform  a  preliminary  inversion  of  real  data.  [1  ]  To  help  interpret 
the  results,  in  parallel  we  invert  data  generated  synthetically  for  a  realistic  geometry  and  set  of 
bottom  parameters. 

The  real  data  that  we  invert  was  obtained  by  G.  Frisk,  J.  Doutt,  and  E.  Hays  in  1981. 
The  associated  experimental  geometry  was  described  in  Section  (1.6)  and  is  presented  again  in 
Figure  (VI.1).  We  win  be  uang  the  data  obtained  from  the  lower  receiver  shown  in  this  figure. 
In  Figure  (VI.2)  we  present  a  velocity  profile  and  density  parameters  for  a  bottom  that  we 
believe  is  comparable  to  the  bottom  where  the  real  data  was  taken.  We  use  this  geometry,  velo¬ 
city  profile,  and  these  density  parameters  to  generate  the  synthetic  data  of  this  chapter.  This 
synthetic  data  is  generated  using  the  hybrid  procedure  described  in  Chapter  IV  and  the  numeri¬ 
cal  Hankel  transform  that  was  described  in  Section  (HI. 7)  }  [2  ]  The  efficiency  of  this  Hankel 
transform  algorithm  made  it  possible  to  obtain  high  quality  results  over  a  large  range  that  would 
otherwise  not  have  been  practical. 

We  begin  by  generating  the  synthetic  data  for  this  geometry  and  bottom.  We  use  the 
numerical  procedure  described  in  Section  (IV. 2)  to  generate  the  plane  wave  reflection  coeffi¬ 
cient,  r(*,),  as  a  function  of  horizontal  wave  number.  Its  magnitude  and  phase  are  presented  in 
Figures  (VI.3a)  and  (VI.3b).  We  see  that  a  pole  is  present  in  r(Jtr)  beyond  the  water  wave 
number.  This  pole  is  due  to  the  low  speed  channel  just  below  the  water-bottom  interface. 
Because  the  source+receiver  height  is  large,  this  pole  will  contribute  an  insignificant  amount  to 


1)  This  algorithm  wai  impiamaamd  la  Foitna  oa  a  VAX-117780  by  Mike  Weagrovitz. 
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Figure  Vl.l  Experimental  geometry  associated  with  the  real  data 
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Figure  VI. 2  The  velocity  profile  and  density  parameters  for  a  bottom  comparable  to  that  where 
the  real  data  was  taken 
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the  pressure  field  and  need  not  be  removed  as  was  done  in  Section  (IV .3b).  Consequently  we 
can  generate  the  field  with  the  hybrid  procedure  described  in  Section  (IV. 3a).  The  magnitude 
and  residual  phase  of  the  associated  synthetic  field  are  presented  in  Figures  (VI.4a)  and  (VI. 4b) 
In  these  figures  very  little  high  frequency  ripple  is  apparent  even  at  large  ranges,  implying  that 
the  field  is  indeed  adequately  represented  and  not  suffering  from  spatial  aliasing. 

Figures  (VI. 5«)  and  (VI. 5b)  present  the  magnitude  and  residual  phase  of  the  synthetic 
field  after  inclusion  of  the  incident  field.  The  regular  behavior  in  these  plots  suggests  that  the 
magnitude  and  readual  phase  are  good  representations  of  the  total  field.  As  further  confirma¬ 
tion  of  the  validity  of  the  total  synthetic  fields  generated  for  this  example,  we  present  the  output 
of  a  ray  program  that  was  run  for  this  profile  in  Figure  (VI.5c).1,2  The  two  synthetic  fields  axe 
in  good  agreement  except  in  the  region  of  the  caustic,  1500m  <  r  <  2000m ,  where  the  ray 
method  is  known  to  be  inaccurate. 

Figures  (VI.6a)  and  (VI.6b)  present  the  magnitude  and  residual  phase  of  the  real  data 
(which  includes  the  source  field).  In  the  region  beyond  the  first  hundred  meters,  the  magnitude 
and  readual  phase  of  die  real  data  behave  regularly,  which  gives  us  confidence  in  them.  The 
interference  pattern  apparent  in  the  magnitude  is  similar  to  that  of  the  synthetic  data.  The  zeros 
in  the  magnitude  are  well  matched  by  the  the  changes  in  the  residual  phase  for  large  ranges.  The 
first  few  hundred  meters  of  the  residual  phase,  however,  looks  significantly  different  from  the 
residual  phase  of  the  synthetic  field.  In  this  region,  changes  in  the  source-height  have  their 
greatest  effect  on  the  measured  field  because  the  geometry  is  most  significantly  affected  by 
source-height  variation  in  this  region.  We  recall  that  the  residual  phase  is  given  by: 

«(r)  -  $(r)  -  *oVr2+(z-ro)2  (1) 

Hie  large  negative  slope  of  the  readual  phase  for  low  ranges  could  be  due  either  to  an  estimate 

of  t0  which  was  too  large,  in  which  case  the  residual  phase  would  display  a  negative  phase 

everywhere,  or  to  an  estimate  of  (r  — xq)2  which  was  too '  large.  We  believe  that  this 

1)  Bat  with  •  lUfhtfy  dMaraet  so  am  fcaifbt  of  125  attars  rtthtr  than  125  mtttrv 

2)  1  wish  to  thank  jin  Doatt  tar  protldiaf  at*  syahotto  Sold. 
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Flf*r»  VI.4a  Log-magnitude  of  the  aynthetic  pr  enure  field  calculated  uaing  the  realiatic  bottom 
parameters  and  geometry 
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Figure  VI.5 b  Residual  phase  of  the  total  synthetic  pressure  field  after  the  inclusion  of  the  source 
field 
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Figure  VI.5e  Log-magnitude  of  the  total  field  calculated  by  a  ray  method 
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uncharacteristic  behavior  is  due  to  imperfect  knowledge  of  the  source  height  in  that  region.  The 
gentle  negative  slope  of  the  unwrapped  phase  for  the  real  data  for  large  ranges  is  probably  due 
to  a  slight  overestimate  of  the  water  wave-number,  kQ. 

Before  we  attempt  to  invert  the  real  experimental  data  to  estimate  the  depth-dependent 
Green’s  function  two  major  factors  must  be  considered.  First,  the  experimental  data  is  available 
only  over  a  finite  range  and  second,  it  is  available  only  at  discrete  points  which  are  not  spaced 
properly  for  our  processing.  The  first  issue  can  be  resolved  by  referring  to  Chapter  V  where  we 
showed  that  for  the  source-height  and  geometry  used  to  obtain  the  experimental  data,  it  was  only 
necessary  to  know  the  field  out  to  about  3040  meters  to  minimize  the  degradation  due  to  win¬ 
dowing.  The  experimental  data  is  available  to  6000  meters.  We  believe,  therefor,  that  window¬ 
ing  should  not  prevent  its  successful  inversion.  The  second  issue  can  also  be  resolved  by  refer¬ 
ence  to  Chapter  V  where  we  showed  that  by  interpolating  the  magnitude  and  unwrapped  phase 
it  was  often  possible  to  translate  the  pressure  field  data  available  on  one  set  of  ranges  to  another. 
We  will  use  the  procedure  developed  there  to  interpolate  the  experimental  data  onto  the  set  of 
ranges  that  we  require  for  processing  by  the  Hankel  transform.  In  parallel  we  will  process  the 
synthetic  data.  The  processed  synthetic  data  provides  a  useful  measure  of  the  success  of  our  pro¬ 
cessing  because  the  depth -dependent  Green's  function  that  we  obtain  can  be  compared  to  the 
true  depth- dependent  Green’s  function  which  is  known  for  the  synthetic  data,  and  presented  in 
Figures  (VL7a)  and  (VI.7b). 

Figures  (VI.8a)  and  (VL8b)  present  the  magnitude  and  phase  of  the  Green’s  function  cal¬ 
culated  by  procesring  die  synthetic  data.  The  synthetic  data  was  originally  available  on  the  grid 
mr  »« •  0, 1,  2,  •  ■  •  .  It  was  linearly  interpolated  (through  its  m»gnitnH»  and  unwrapped 

phase  as  described  in  Section  (V.3e))  onto  the  grid  required  for  processing,  —  for 

A 

n  »  0, 1,  2,  •  •  •  with  A  m  1.2  and  where  K,,  n  ■  0, 1,  2,  •  •  •  are  die  zeros  of /<j(x). 

The  agreement  between  the  estimate  of  the  synthetic  Green’s  function  obtained  by  process¬ 
ing  die  synthetic  field  and  shown  in  Figures  (VI. 8a)  and  (VI.8b)  and  the  true  Green’s  function 


Figure  VI.7a 


t 


Figure  VI. 7b  Ruse  of  the  true  depth-dependent  Green's  function  for  the  synthetic  data 
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for  the  synthetic  data  shown  in  Figures  (VI.7a)  and  (VL7b)  is  excellent,  particularly  the  agree¬ 
ment  in  the  magnitudes.  The  phases  differ  slightly  for  low  horizontal  wave  numbers.  We  believe 
that  this  is  due  to  small  errors  in  the  synthetic  field  for  low  ranges.  The  phases  differ  dramati¬ 
cally  in  the  evanescent  region  beyond  the  water  wave  number,  where  the  magnitude  of  the 
Green’s  function  is  very  small  and  consequently  the  phase  is  probably  dominated  by  noise.  The 
agreement  in  general  between  the  true  Green’s  function  for  the  synthetic  data  and  the  Green’s 
function  estimated  from  die  synthetic  data  is  excellent,  however,  and  confirms  the  results  of 
Chapter  V  which  indicated  that  for  the  sampling  rate  and  range  of  values  over  which  the  data  is 
known,  it  should  be  possible  to  determine  the  depth-dependent  Green’s  function. 


Figures  (VL9a)  and  (VL9b)  present  the  magnitude  and  phase  of  die  Green’s  function  cal¬ 
culated  from  the  real  data.  Except  for  low  wave  numbers,  the  magnitude  of  this  Green’s  function 
displays  many  of  the  features  of  the  synthetic  Green’s  function,  including  die  same  overall 


envelope  due  to  die  — 7— » 

v» Fj? 


source  spectrum  term,  and  the  interference  pattern  arising  from 


the  interaction  of  that  portion  of  the  Green's  function  associated  with  the  source  and  that  portion 
associated  with  die  reflected  field.  The  total  Green’s  function  also  decays  rapidly  at  the  water 

wave  number,  as  it  should  due  the  the  migration  terms.  In  the  evanescent 

region,  kr>  ig,  we  see  only  noise,  comparable  to  the  noise  we  see  superimposed  upon  the  rest 
of  the  spectrum. 

At  low  horizontal  wave  numbers  the  Green’s  function  for  the  real  data  does  not  look  like 
die  Green’s  function  for  die  synthetic  data.  Very  near  die  origin  we  see  a  large  peak  not 
apparent  in  the  total  Green’s  function  for  the  synthetic  data.  This  peak  is  probably  due  to  con- 
centrstion  of  noise  power  there  by  the  Hankei  transform  as  discussed  in  Section  (Q.8).  For 

sightly  larger  wave  numbers  the  magnitude  displays  a  jagged  appearance  not  seen  in  the  total 

% 

Or  sen’s  function  for  the  synthetic  data.  In  this  region,  the  stationary  phase  approximation  for 
the  Sommerfeid  integral  is  fairly  good,  allowing  us  to  associate  the  behavior  of  the  Green’s  func¬ 
tion  at  low  horizontal  wave  numbers  with  die  behavior  of  the  pressure  field  at  tow  ranges.  The 
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Figure  VI.9a  Phase  of  the  depth-dependent  Green’s  function  estimated  by  talcing  the  Henkel 
transform  of  the  interpolated  real  data 
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uncharacteristic  behavior  of  the  Green's  function  at  low  horizontal  wave  numbers  is  consistent 
with  the  uncharacteristic  behavior  of  the  residual  phase  that  we  observed  for  low  ranges  and  may 
be  doe  to  variations  in  the  source-height.  Some  of  this  apparent  jitter  in  the  Green’s  function  of 
the  real  data  may  be  due  in  part  to  variation  in  the  source-height.  A  rough  sampling  of  die 
source-height  over  the  course  of  the  experiment  was  available  from  die  experimental  records. 
We  interpolated  between  available  amples  using  splines  to  obtain  a  rough  estimate  of  die 
source-height  variation  present  during  the  course  of  die  experiment  The  result  is  presented  in 
Figure  (VL10).  This  curve  is  sufficiently  similar  to  the  sum  of  the  two  low  frequency  cosines  dis¬ 
cussed  in  Section  (V.5b.iii)  to  qualitatively  interpret  die  effect  of  source-height  variation  for  this 
experiment  in  term  of  the  results  presented  there.  The  analysis  of  Section  (V.5)  shows  that 
sinusoidal  variation  in  die  source-height  causes  the  estimated  Green's  function  to  be  a  rever¬ 
berant  version  of  the  true  Green's  function,  particularly  for  low  kr  corresponding  to  large  kg . 
Because  the  frequency  of  the  variation  is  very  small,  the  main  effect  is  to  anear  the  ewimate  of 
y/krG  (kr).  As  stated  in  that  section,  the  phase  of  the  estimated  Green's  function  might  be  more 
seriously  corrupted  than  its  magnitude.  The  phase  of  the  depth-dependent  Green's  function 
estimated  from  the  real  data  and  shown  in  Figure  (VI.  9b)  does  not  strongly  resemble  the  phase 
of  the  synthetic  Green's  function.  The  overall  good  appearance  of  the  magnitude  of  the  total 
Green's  function  and  dm  poor  appearance  of  it  phase  is  consistent  with  the  the  degradation  that 
would  be  expected  from  source-height  variation. 

Figures  (VLlla)  and  (VI. lib)  riiow  the  magnitude  and  phase  of  the  plane  wave  reflection 
coefficient  generated  from  the  Green's  function  calculated  from  the  synthetic  data  and  shown  in 
Figures  (VL8a)  and  (VI.8b).  Figures  (VI.12a)  and  (VI.  12b)  present  the  magnitude  and  phase 
of  the  plane  wave  reflection  coefficient  calculated  from  the  Green's  function  for  the  real  date. 
The  estimate  for  the  reflection  coefficient  for  the  real  data  does  not  appear  to  be  a  good  one  at 


Because  the  plane  wave  reflection  coefficient  is  obtained  from  the  total  depth -dependent 
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Figure  VI. lie  Magnitude  of  the  plane  wave  reflection  coefficient  generated  from  the  depth- 
dependent  Green's  function  estimated  from  the  synthetic  data 
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Figure  VI. 12a  Magnitude  of  the  plane  wave  reflection  coefficient  generated  from  die  depth- 
dependent  Green's  function  estimated  from  the  real  data 
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Figare  VI. 12b  Phase  of  the  plane  wave  reflection  coefficient  generated 
dependent  Green’s  function  estimated  from  the  synthetic  data 
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Green’s  function  by  first  coherently  subtracting  the  source  contribution  and  then  muidpiying  by  a 

term  with  a  rapidly  varying  phase  (e  ~k'  '* _r#l  )  errors  in  the  phase  of  the  total  Green’s 
function  would  seriously  degrade  the  estimate  of  plane  wave  reflection  coefficient.  The  estimate 
for  the  reflection  coefficient  is  probably  much  worse  the  estimate  for  the  Green’s  function 
because  of  the  phase  errors  in  the  estimate  for  die  total  depth-dependent  Green’s  function. 

In  canduston,  we  believe  most  of  the  error  apparent  in  the  Green's  function  for  die  real 
data  to  be  due  to  variation  in  the  source-height  near  die  origin.  Direct  evidence  of  this  is  the 
anomalous  residual  phase  variation  in  the  region  r<  300  meters.  The  error  in  the  estimated 
Green’s  function  for  very  small  horizontal  wave  numbers  is  probably  due  to  additive  noise.  The 
errors  in  the  estimator!  reflection  coefficient  are  probably  due  to  imperfect  knowledge  of  the 
source-receiver  geometry  that  affects  the  coherent  additions.  Overall,  however,  we  are  greatly 
encouraged  by  die  good  appearance  of  die  magnitude  of  the  total  depth-dependent  Green’s  func¬ 
tion  determined  from  the  real  data.  The  interference  structure  and  the  overall  envelope  suggest 
that  we  are  very  dose  to  being  able  to  estimate  the  plane  wave  reflection  coeffident  from  real 
data.  Work  still  needs  to  be  done  to  compensate  for  the  effect  of  source-height  variation. 

The  potential  returns  from  the  successful  inversion  of  pressure  field  data  to  obtain  the  plane 
wave  reflection  coeffident  are  enormous.  Such  a  successful  inversion  is  a  vital  step  in  the  process 
of  inferring  the  physical  parameters  of  the  bottom  from  acoustic  measurements.  [3, 4  ]  The  abil¬ 
ity  to  make  such  inferences  is  of  great  interest  to  oceanographers  and  to  exploration  geophysi¬ 
cists.  A  successful  inversion  would  also  make  it  posable  to  predict  the  fields  assodated  with  an 
arbitrary  source-receiver  geometry  from  one  set  of  measurements.  This  would  greatly  facilitate 
acoustic  imaging  in  the  ocean. 
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CHAPTER  VII: 

CONTRIBUTIONS  AND  FUTURE  WORK 

Vn.l)  Contributions 

In  Hus  thesis  we  have  studied  both  the  numerical  generation  of  synthetic  pressure  fields 
from  dm  plane  wave  reflection  coefficient  and  the  inversion  of  measured  prewure  field  data  to 
estimate  die  plane  wave  reflection  coefficient  We  developed  and  implemented  algorithms  that 
efficiently  generate  high  quality  synthetic  fields.  We  studied  the  major  issues  affecting  the  inver¬ 
sion  of  experimental  data  and  were  able  to  the  depth-dependent  Green's  function  from 

measured  data  taken  in  the  ocean  with  a  high  degree  of  success.  We  isolated  source-height  varia¬ 
tion  as  a  major  factor  preventing  the  successful  estimation  of  the  plane  wave  reflection  coefficient 
it  this  time. 

As  a  foundation  for  our  studies  we  explored  the  Hank  el  transform  in  depth.  In  Chapter  II 
we  derived  a  number  of  important  properties  including  the  effects  that  windowing  and  sampling 
a  function  have  on  its  Hankel  transform.  Our  sampling  results  show  that  the  associated  degrada¬ 
tion  is  often  a  more  severe  problem  for  die  Hankel  transform  than  for  the  Fourier  transform.  In 

particular  it  can  seriously  degrade  synthetically  generated  pressure  fields  which  decay  as  or 
even  more  slowly  and  its  effect  should  always  be  carefully  considered. 

In  Chapter  H  we  also  studied  dm  noise  properties  of  the  Hankel  transform.  We  showed 
that  if  a  function  is  sampled  on  a  square  root  grid  in  a  noisy  environment,  its  Hankel  transform 
will  have  superior  noise  properties  more  characteristic  of  the  underlying  two  dimension  Fourier 
transform  which  the  Hankel  transform  represents  in  dm  presence  of  cylindrical  symmetry. 

In  Chapter  IH  we  considered  a  number  of  numerical  techniques  for  performing  the  Hankel 
transform.  We  presented  new  results  strengthening  existing  procedures  such  as  dm  asymptotic 
and  backanear  methods  as  well  as  an  efficient,  exact  method  developed  as  part  of  this  diesis. 

hi  our  development  of  algorithms  to  generate  high  quality  synthetic  date  we  presented  a 
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number  of  hybrid  numerical-analytical  techniques  that  (ready  improve  the  quality  of  synthetic 
data.  In  the  course  of  developing  a  technique  that  can  adequately  handle  the  effects  of  guided 
modes  in  slow  speed  layers  under  the  ocean  bottom,  we  derived  an  expression  that  may  be  of  use 
for  developing  field  expressions  in  modal  expansions  required  to  be  accurate  in  both  the  near 
and  far  fields.  We  also  presented  a  well  behaved  numerical  procedure  for  implementing  the 
Thomson-HaskeU  approach  for  generating  die  plane  wave  reflection  coefficient. 

hi  Chapter  V  we  developed  the  major  issues  affecting  die  inversion  of  measured  field  to 
obtain  the  plane  wave  reflection  coefficient  On  the  baas  of  this  development  we  were  able  to 
identify  the  sources  of  error  in  an  actual  inversion.  The  phase  unwrapping  and  interpolation 
results  presented  in  this  chapter  also  significantly  improved  the  results  of  die  processing  of  the 
experimental  data  in  Chapter  VI. 

In  Chapter  VI  we  performed  a  preliminary  inversion  of  real  data  to  obtain  estimates  for 
die  depth-dependent  Green’s  function  and  the  plane  wave  reflection  coefficient  The  results 
presented  in  this  chapter  represent  a  significant  advance  towards  the  complete  inversion  of  meas¬ 
ured  pressure  field  data  to  obtain  the  plane  wave  reflection  coefficient.  We  were  able  to  generate 
a  good  estimate  for  the  depth-dependent  Green’s  function  and  were  able  to  associate  the  effects 
of  source-height  variation  with  the  degradation  in  the  estimate  for  the  plane  wave  reflection 
coefficient. 

At  this  point  work  is  continuing  towards  the  complete  estimation  of  the  plane  wave  reflec¬ 
tion  coefficient  from  real  data.  The  foundations  laid  by  the  work  presented  in  this  thesis  provide 
a  strong  base  for  future  work  in  this  area.  In  addition  they  suggest  research  in  a  number  of 
related  areas.  Some  of  thew  are  prewnted  in  the  next  section. 
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VH.2)  Future  work 


i)  Cylindrical  to  Cartesian  Coordinate  Systems 


la  this  thesis  we  have  dealt  with  problems  cast  ia  a  cylindrical  coordinate  system.  Ia  that 


eoordhute  system  the  familiar  Former  transform  of 


systems 


die  leas  familiar 


Hankel  trmasfbrm.  Ia.  that  coordinate  system  the  deaa  properties  of  additive  white  Gaussian 
noise  through  the  Fourier  transform  were  obscured  obscured  until  a  square  root  grid  was  intro¬ 
duced.  In  that  coordinate  system  windowing  and  aliasing  approximately  affected  vr/(r) 

instead  of  f(r).  In  that  coordinate  system  the  familiar  impulse  8(jr)  became  ^  .  In  that 

a 2  ,  a2  I 

coordinate  system  the  — —  operator  became  V2  m  — -  +  — —  so  that  the  operator  which 
Hx2  Hr2  r  9r 


nulls  a  pole  in  a  cartesian  coordinate  system: 


m  \ih  +  m  "  ~*(x)  (1) 


p. 


became  the  less  familiar: 


[v2  +  P,2l/(r)  -  [v2  +  p 2 1 J — o(Pr )p^  P  *  — ^-Jo(pr)pdp  -  -- (2) 

L  r  L  lorPi  0  p  -P/  r 

In  diort  we  frequently  found  that  familiar  problems  in  a  cartesian  system  became  more  difficult 
when  caste  in  a  cylindrical  coordinate  system.  The  reverse  is  also  true,  however.  In  Section 
(m.7)  we  developed  an  efficient  numerical  algorithm  for  the  Hankel  transform  by  mapping  it 
into  a  Fourier  transform.  The  mapping  was  accomplished  with  the  Abel  transform: 
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In  Section  (HI. 7)  we  developed  an  efficient  numerical  algorithm  for  evaluating  the  Abel 
tranafonn.  One  exciting  area  of  future  research  is  the  extend  on  of  maximum  entropy  and  other 
spectral  estimation  techniques  into  die  cylindrical  domain  as  we  now  describe. 

An  estimation  scheme  that  might  be  of  value  for  estimating  the  plane  wave  reflection  coef¬ 
ficient  is  illustrated  in  Figure  (VH.2a.l).  Instead  of  estimating  the  plane  wave  reflection  coeffi¬ 
cient  directly  we  estimate  the  portion  and  residue  of  its  poles  (and  possibly  aeroa)  in  the  complex 
plane.  We  do  this  because  the  •wtmarinw  of  parameters  instead  of  a  function  is  a  much  better 
posed  problem  when  the  agnal  available  for  analysis  has  been  corrupted.  Because  the  plane  wave 
reflection  coefficient  is  related  to  the  measured  field  by  the  (cylindrical)  Hankel  transform  and 
not  the  Fourier  transform,  spectral  estimation  techniques  available  in  the  literature  of  digital  sig¬ 
nal  processing  do  not  apply.  If  we  first  process  the  pr enure  field  data  with  an  Abel  transform, 
however,  the  resulting  signal  has  die  same  poles  and  zeros  but  now  in  its  Fourier  transform. 
Modern  spectral  estimation  techniques  can  therefore  be  used.  We  introduce  the  caution  that  the 
effect  that  branch  cuts  have  on  this  procedure  must  be  studied  with  care. 

b)  Analytical-Numerical  Algorithms 

The  hybrid  analytical-numerical  technique  used  to  implement  the  Abel  transform  in  Section 
(IH.6)  is  a  vary  general  procedure  and  springs  from  classical  numerical  methods  of  long  standing. 
Traditionally,  difficult  integrals  are  evaluated  numerically  by  removing  their  singular  behavior  as 
much  as  posable  through  coordinate  changes  and  change*  of  variables  and  then  numerically 
transforming  die  result.  The  success  of  the  hybrid  method  points  out  that  in  fact  it  is  often  desir¬ 
able  to  do  just  dm  opposite.  The  integral  should  be  manipulated  to  produce  as  much  angular 
behavior  as  possible.  The  singularities  can  be  integrated  analytically  and  will  not  suffer  from 
numerical  degradation.  If  dm  angularities  are  removed  properly,  dm  remaining  numerical  por¬ 
tion  of  dm  integral  will  be  well  behaved  where  it  dominates  and  subordinate  to  the  analytically 
determined  portions  of  the  integral  where  it  does  not  Hie  art  in  this  procedure  is  casting  the 
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MEM 

SPECTRAL 

ESTIMATION 


Figure  VII.2a.l  The  Proposed  method  silo  wing  MEM  to  be  used  for  problems  caste  in  cylindri¬ 
cal  coordinates 
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an&lyticaHy  determined  portions  of  the  integral  carefully  to  insure  that  die  numerical  portion 
does  not  have  infinities  or  undesirable  asymptotics  to  cancel.  Had  the  e~br  factor  not  been 
included  in  the  procedure  of  Section  (HI.  6),  for  example,  the  hybrid  method  for  the  Abel 
transform  would  not  have  worked. 

Manipulating  functions  so  that  they  can  best  be  represented  by  parameterized  functions  and 
samples  is  part  of  die  general  issue  of  computer  representation.  As  software  systems  become 
smarter  this  kind  of  approach  will  become  increasingly  more  important. 

c)  Waveguides 

In  the  course  of  generating  synthetic  data,  we  evaluated  the  integral: 

*  . 

Hr* ,Pi)  *  f "7™-' -r i^^1* '/oCpOp* p  (i) 

0  p^-pf 

and  showed  that  it  satisfied 

[Vr+P*]/(»-.*,P<)  “  *  0  (2) 

We  associated  I(r  ,z  ,p()  with  the  contribution  of  the  pole  at  p,  because  our  integrals 

always  included  the  terms  -■>  * . I  M  well.  The  advantage  of  this  formulation  was 

V*2_p2 

that  l(r  ,z  ,p()  is  everywhere  finite,  even  at  r  =  0.  The  classical  contribution  associated  with  a 
pole  is 

“Y^l)(p/0  when  Im(P/)>0  (3) 

The  Hank  el  function  above  has  a  logarithmic  singularity  as  the  origin.  Physically,  poles  in 

the  depth-dependent  Green’s  function  make  only  finite  contributions  to  die  field.  The  migration 

term  *<v*J-p*I*  I  windows  die  pole  in  the  Hankel  domain  so  that  its  contribution  to  the  pressure 

field  is  everywhere  finite.  For  this  reason  we  included  the  migration  term  into  our  pole  expres- 

sion.  For  convenience  we  also  included  the  source  term  * 

V*2_p2 
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Our  formulation  was  carefully  constructed  to  insure  that  the  numerics  were  not  required  to 
generate  infinities.  It  is  potentially  useful  for  many  problems  other  than  those  amndered 
directly  in  this  thesis.  One  set  of  problems  concerns  the  calculation  of  fields  inside  (pomibiy 
leaky)  waveguides. 

We  develop  an  espremion  for  die  field  inside  a  (dielectric)  waveguide  ariung  from  a  point 
source  using  die  plane  wave  formulation  of  this  theas.1  Figure  (VII. 2c.  1)  presents  the  geometry 
of  the  waveguide  and  the  waves  present 

The  radial  and  time  variation  of  all  fields  is  given  by: 


J^rr)e~^  (4) 

and  will  be  suppressed.  We  will  use  $  for  the  vertical  wave  number,  and  kr  are  related 

through 


k2  -  P2+i,2  (5) 

The  source  field  is  given  by  and  is  the  portion  of  the  field  that  would  be 

present  even  in  die  absence  of  impedance  contrasts.  We  specify  the  boundary  conditions  at 
z  *  h  by  giving  the  plane  wave  reflection  coefficient  there,  r7(kr),  and  at  z  =  0  by  TB(kr). 
These  two  interfaces  together  give  rise  to  an  up-going  wave,  P+e'P*.  and  a  down-going  wave 
that  would  not  be  present  without  the  impedance  contrasts  in  the  regions  z  a h  and 
z  SO.  The  total  up-going  and  down-going  fields  for  zq<z  jS h  is  given  by: 


P+e'**  +f»/e'W,",a)  UP 
P-e'1*1  DOWN 

The  plane  wave  reflection  coefficient  at  z  “h  provides  the  boundary  condition: 

r 

In  the  region  OSz  <z0  we  have 

P  -e  +  P/e'Wr,_,)  DOWN 
P+***  UP 


(6) 

CO 

(8) 


A  dwissdoa  tor  mo  field  fauide  web  a  waveguide  an  also  bo  toaad  is  [1 ) 
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and 


P+elfi0 


P  + 


P-e’W  +  Pte P.  +  P,e‘*i0 


(9) 


Pj  is  die  known  point  source  strength 


=  Tt  and  TB  are  given  as  wefl. 


The  two  boundary  conditions  are  sufficient  to  determine  P  +  and  P  _  from  these  quantities.  We 
write: 


rrP+e1**  +  -P-e^ 

r bp.  +  rBpte iUt  ~p+ 

Solving  for  P  +  and  P  -  we  have 


(10) 


and 


«-'«*-*•>  -  rreip(*-Ifl)k 

P+  *  J 

-f-e-'P*  -I>*'P* 
rs 

Ir,*1™**  -  e,p(*-,o)k 

/»_  =  -L - - - J —  *  C  JPj 

-T.e'P* 


(ID 


(12) 


The  total  field  in  the  waveguide  (  0<z  <h  )  is  given  by 


P(r,z)  *  / 
o 

with 


[e'V*J"*'Il«'»0l  +  c+r'V*,-*,J»  +  C  e-i V*2-*,*. 


(13) 


C+ 


-i-e-'P*  -  rre'P* 

r* 


C_ 


|^0(*+»o)  _  £<P(*-*o) 


(14) 


(15) 


-rte+,p* 

*r 

The  zeros  of  e  _,P*  -  Tr^P*  -  0  contribute  poles  to  die  depth -dependent  Green's 
function  and  give  rise  to  die  modes  of  the  waveguide.  Each  of  these  poles  in  Equation  (13) 
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makes  a  contribution  to  the  field  of  the  fora 

K 

I(r,z,pt)  =  J-4— T€|V*5=^l'l-3^-T/0(pr)pdp  (16) 

0  V*2-p2  p2-pr 

and  the  development  we  used  for  calculating  fields  in  the  presence  of  poles  applies. 

It  should  be  noted  that  the  inverse  problem,  that  of  resolving  modes  in  a  waveguide,  can  be 
cast  into  die  damical  signal  processing  fora  of  finding  poles  in  die  Fourier  transform  by  first 
generating  the  Abel  transform  of  the  pressure  field!  The  effect  of  brands  Hb«  on  this  approach 
needs  to  be  studied,  however. 

It  is  also  possible  to  construct  nulling  operators  to  estimate  die  pole  positions  (as  is  done  in 
Maximum  Entropy  spectral  estimation  techniques)  by  using  Equation  (2),  which  does  indude 
some  branch  line  effects. 
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(8) 


jto  S22^co.^L  -  2(-D‘5(f  “  2*) 

aniT-^-  * 

2 

This  is  the  result  needed  in  the  text. 

To  evaluate  the  limit  of  Equation  (2): 

cotNvx  .  vx 

4 

we  first  consider  die  behavior  of  near  x  —  0.  The  behavior  at  the  other  zeros  of  sin-^- 

**— 


(9) 


will  be  similar.  We  assume  that  x  is  sufficiently  small  so  that  sm-~p-  can  be  replaced  the  first 


term  of  its  Taylor  series,  We  consider: 


limSSfitZl 

*-»  *X 
2 


(10) 


Instead  of  evaluating  the  limit  direcdy,  we  look  at  the  Fourier  transform  of  the  limit.  We 


evaluate  this  by  taking  die  limit  of  the  Fourier  transforms  of  each  term  and  write:1 


lim 

W-* 

cosN  vx 

-lim  FT 

coaNvx 

It* 

vx 

k 

2 

* 

2 

l 

(11) 


The  Fourier  transform  of  c0^,rx  can  be  found  by  convolving  the  Fourier  transform  of 

T 


cos N  vx  with  the  Fourier  transform  of  — . 

vx 


so  that: 


FT 


etnNvx 


vx 

2 


(12) 


lynasa  atpt  am  b«  rigorously  jostifiad  by  arise  gaaar «Mart  Auettoa*.  [1  ]  Thay  prataaa  that  two  fa aw 
tloas  an  said  to  ba  aqual  if  tha  mult  of  oonroMng  thafr  dMt—ca  with  toy  bad-Usriwd  fuacdoa  is  ah 
ways  san.  Akamataiy,  two  faacdoas  an  irid  to  ba  aqual  if  tha  Fouriar  transform  of  thaii  —  is 
swo  lor  my  flate  bead. 
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APPENDIX  II: 

THE  VALUE  OF  THE  KERNEL 

FOR  THE  NUMERICAL  PORTION  OF  THE  HYBRID  ALGORITHM 
AT  THE  WATER  WAVE  NUMBER 

la  this  appendix  we  derive  the  value  of  £(fto)  discussed  la  Section  (IV.3e).  L(k,)  is 


defined  by  Equation  (IV.3a.15)  as: 


*<*,) 

and  we  seek  to  evaluate  the  limit 


(1) 


lim  L(kr) 

*r~*  0 

under  the  condition  that  the  impedance  of  the  bottom,  Z(kr),  is  finite  at  kt  =  k0. 


(2) 


At  kr  =  Jt o  Equation  (1)  takes  the  indeterminate  form  — .  We  evaluate  die  limit  (2)  by 
using  L'  Hospital’s  rule: 

i  |r(*,)  -  1  -±-t  |r(*,)  -  r(*0)  j*' V^l‘ 1 

w ' 


lim 


Km 

*,-*o 


dlr  Vi*  ~k' 


(3) 


After  separating  out  the  terms  that  approach  zero  as  this  expression  becomes: 


Um 


L(kr)  -  lim  ~-'2  r(kr)e,Vk*~i'li‘ 1 


Mi 


(4) 


We  now  express  t(kr)  in  terms  of  the  characteristic  impedance  of  the  upper  half  space,  Zfc 
and  the  impedance  at  the  interface  which  we  will  denote  at  Z\.  Both  Zq  and  Z}  are  of 

k r  in  general.  In  terms  of  them  r(*,)  is  given  by: 

Zi-Zo 


Taking  derivatives  we  find: 


IX*,) 


t(kr)  -  2 


z,+z, 


Z)Zc  ~  2tZ| 


(Z,+Zo)a 


(5) 


(«) 
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We  now  use  the  characteristic  impedance  of  the  upper  half  space: 

Z0 


VkFtf 


(7) 


where  pg  is  die  density  of  the  half  space,  u  is  the  temporal  frequency,  and  k0  m  —  where  Co  is 

Co 


die  qjeed  of  sound  in  the  typer  half  space.  Substituting  Equation  (7)  into  Equation  (6)  and 
evaluating  ZQ  we  have: 


f(*r)-2 


i  "0-  _  »rW  J 

T  n  ,  f'"8**  T 

_  4 

**"  ‘ 

*  i 

zUki-kfi+iz^Vki-k?  *  f 

?  ViJ-t,*  ki~k? 

(«) 


Substituting  (8)  into  (4)  we  find 


„  v  ..  -iVkf-*}  , 
lim  £(*,)  ■  lim - - - 


.i  *rPO<*>  _ 


Zhki  “k,?)+2Z  iPo®  V ki  -kr2+p$v>2 


(9) 


or 


L(k0) 


-2JZt(k0) 


(10) 


provided  that  lim  Vk^  -  k^Z  i(k, )  *  0 

*,-*o 


If  the  interface  is  between  two  isovelocity  half  spaces,  the  expression  (10)  for  L  (*0)  on  be 

pl(0 

written  directly  in  terms  of  die  material  parameters.  For  this  case  Zt(kr)  *  ^  i '  z,(*°)  “ 


finite  because  k j  #fcg  ( if  kj  *  kg  there  would  be  no  interface).  L  (kg)  i*  given  by 
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APPENDIX  ni: 

EVALUATION  OF  THE  POLE  CONTRIBUTION  TO  THE  FIELD  FOR  SECTION  (IV.3B) 
Here  we  evaluate  the  pole  contribution  to  Equation  (1): 


/(»*,*;*)  -  jf-r-TTTT- 7  e'V*3I?|,|/o(pOprfp  Cl) 

0  r-pi  V*2_p2 

We  evaluate  Equation  (1)  by  determining  a  partial  differential  equation  that  it  satisfies  and  solv¬ 
ing  that  equation. 

Taking  die  second  partial  derivative  of  Equation  (1)  with  respect  to  z  we  have 


j£/ -  j£g- 


1'-*Jfrr)pdpW)S -Ji-tf- 


'■/«(*•  )p*p  (2) 


If  we  use 


8(z)/(z)  =  8(z)/(0)  for  any  /(r) 


then  Equation  (2)  becomes 


rrHr s #,)  -  jr^z4r7n-T''v*i"pJ|‘  (*) 

0  p  -P(  v*2- p2  o  P  -Pi 

Putting  it  all  together  we  have: 

-  (p*  -  k 2)  I(r,z;pi)  -  f  \Zk2-p2  l/«Cp»’)pdp  ”  2&(r)  J^2  *  piJo(pOpd(i*) 

If  we  define  02  *  p2  —  t2  ,  choosing  real  part  of  P  >  0,  and  use 


* 

/  a1  i^oCpr^P  -  *•(“*«»  ^"(Pi)!  ipir)  -  —p-BfPfar)  when  /»(p,)  >0  (I) 

•  P  “Pi  * 

then  (A2.5)  becomes 

|^,-P2J/(»‘^y><)  -  ~r2+g2"  +  f»8(r)ffrf1)(p|r)  (I) 

when  /m(p()  >  0. 
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The  Green’s  function  for  this  differential  equation  is  given  by 

C(r.,.i)  -  W 

Using  this  for  the  impulsive  response  and  convolving  with  the  continuous  driving  function  to 
rthtain  the  particular  solution  we  obtain: 

When  An  (pi)  >  0. 

The  general  form  of  this  expreaion  which  is  valid  for  all  p<  is  given  by: 


(11) 
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